A protein-ligand interaction prediction method and optimization method

Abstract

The present application belongs to the technical field of computer-aided drug design, and discloses a protein-ligand interaction prediction method and an optimization method. The method comprises the following steps: constructing a data set and making label data; establishing a protein-ligand complex model, encoding the protein-ligand into a protein binding pocket graph and a protein-ligand atom interaction graph, and the protein-ligand atom interaction graph is an isomorphism graph; the isomorphism graph comprises two types of nodes and four types of edges; constructing a scoring model based on an edge-based graph attention network, comprising a feature extraction module and two decoders; training the scoring model; and using the trained scoring model to predict the ligand binding pose root mean square deviation RMSD and the protein-ligand binding strength pkd. The scoring model can simultaneously predict the ligand binding pose root mean square deviation RMSD and the binding strength pkd, and can simultaneously have excellent performance in scoring, ranking, docking and screening tasks, thereby improving the universality and transferability of the model.

Description

Technical Field

[0001] This invention belongs to the field of computer-aided drug design technology, specifically relating to a protein-ligand interaction prediction and optimization method, and more particularly to a protein-ligand interaction prediction and optimization method based on geometric graph neural networks. Background Technology

[0002] A key task in computer-aided drug design is the discovery of lead compounds with high binding affinity to pathogenic proteins. Molecular docking and virtual screening are crucial technologies in computer-aided drug design, with molecular docking being a core technology of virtual screening. Since the prediction of binding affinity largely depends on the position of the ligand in the binding pocket, selecting the optimal docking posture from numerous options is essential, and this also affects the reliability of the virtual screening results.

[0003] Docking is the process of attaching an unfamiliar ligand molecule to the binding pocket of a protein, resulting in the specific ligand conformation (docking pose) of the ligand molecule within the binding pocket. Docking can be divided into two steps: sampling and scoring.

[0004] (1) Sampling: This refers to searching for the specific location and orientation of small molecule binding within the protein pocket. Current docking software typically uses random sampling algorithms (such as the Genetic Algorithm and Monte Carlo algorithms) to search for the binding posture of ligands within the protein binding pocket;

[0005] (2) Scoring: This refers to predicting the binding energy between proteins and ligands in different binding postures using a scoring function. Most current molecular docking software uses traditional scoring functions, such as physics-based, empirical, and knowledge-based functions. However, traditional scoring functions are not ideal in predicting protein-ligand interactions. In recent years, applying machine learning (ML) and deep learning (DL) algorithms to predict protein-ligand interactions has become a trend. Table 1 below shows commonly used scoring functions, which have been applied to re-scoring and ranking molecular binding postures, such as in GNINA.

[0006] Table 1. Commonly Used Scoring Functions

[0007] Model Name Scoring function type RTMScore Deep learning DeepRMSD Deep learning AutoDock Vina Experience-based scoring function <![CDATA[G NINA 1.0]]> Deep learning

[0008] Current models for assessing drug binding posture and strength are still limited by accuracy and generalization ability, resulting in unsatisfactory performance of traditional scoring functions in predicting protein-ligand interactions. Specifically, most scoring functions do not provide true physical meaning in their predictions, and they are typically only applicable to specific scenarios, making them difficult to transfer to other applications. For example, Pafnucy can predict the binding affinity of protein-ligand crystal structures, but performs poorly in docking and screening tasks; RTMScore performs excellently in docking and screening tasks, but performs poorly in scoring and ranking tasks, and its absolute scores have no real physical meaning. Currently, models that are applicable to scoring, ranking, docking, and screening tasks simultaneously are almost non-existent.

[0009] Another major challenge in current virtual screening is the flexible docking of ligands. When docking an unfamiliar small molecule into the protein target pocket, current molecular docking programs struggle to obtain a precise ligand conformation (docking pose). The goal of molecular docking is to generate a ligand binding posture within the protein binding pocket that closely approximates the native conformation. A commonly used metric for evaluating the quality of binding postures generated by docking software is the root mean square deviation (RMSD). A larger RMSD indicates a greater difference between the binding posture and the native conformation, and vice versa.

[0010] There are two main reasons why the accuracy of molecular docking is currently limited:

[0011] (1) Current molecular docking assumes that proteins are rigid bodies and only considers the flexibility of small molecules. However, in the actual process of small molecules binding to proteins, both molecules will change their structures simultaneously due to interaction and mutual influence until they reach the most stable conformation. Therefore, this fundamentally makes the structure itself less realistic.

[0012] (2) The scoring function is not precise enough to capture key interaction information and subtle differences in the structure. Summary of the Invention

[0013] This invention aims to solve at least one technical problem existing in the background art, and provides a protein-ligand interaction prediction method and optimization method, so that the scoring model can simultaneously predict the root mean square deviation of ligand binding posture and the binding strength pkd, and can have excellent performance in scoring, ranking, docking and screening tasks, thereby improving the universality and transferability of the model in this field.

[0014] To achieve the above technical objectives, the present invention adopts the following technical solution:

[0015] A method for predicting protein-ligand interactions, the method comprising the following steps:

[0016] Step S1: Construct the dataset and create labeled data;

[0017] Step S2: Establish a protein-ligand complex model;

[0018] The protein-ligand relationship is encoded into two diagrams: a protein binding pocket diagram and a protein-ligand atom interaction diagram.

[0019] The protein-ligand atomic interaction diagram is an isomer diagram of protein-ligand interactions, G. RR,LL,RL The heterogeneous graph G RR,LL,RL There are two types of nodes and four types of edges in the system;

[0020] The two types of nodes are: protein atom nodes. and ligand atomic nodes

[0021] The four types of edges are: white matter node - protein node. Ligand node - ligand node Protein node-ligand node Ligand nodes - protein nodes

[0022] Step S3: Construct an edge-based graph attention network scoring model;

[0023] The scoring model includes: a feature extraction module, an RMSD decoder, and a pkd decoder;

[0024] The feature extraction module includes two independent edge-based graph attention layers for encoding features from the protein binding pocket graph and the protein-ligand atom interaction graph into a latent space for protein-ligand interactions.

[0025] The RMSD decoder and the pkd decoder are used to decode the latent space to obtain the predicted root mean square deviation of ligand binding attitude (RMSD) and protein-ligand binding strength (pkd), respectively.

[0026] Step S4: Train the scoring model;

[0027] The two graphs encoded in step S2 are input into the scoring model constructed in step S3. After feature extraction, encoding and decoding of the protein-ligand interaction latent space, backpropagation and parameter update, the trained scoring model is obtained.

[0028] Step S5: Using the trained scoring model, predict the interaction between the protein and ligand to be tested, and obtain the predicted root mean square deviation of ligand binding posture (RMSD) and protein-ligand binding strength (pkd).

[0029] Furthermore, the dataset includes a validation set, a test set, and a training set, and generates more virtual binding postures for the ligands in the training and validation sets through molecular docking, thereby expanding the sample data of the training and validation sets.

[0030] The dataset includes RMSD label data and pkd label data;

[0031] The root mean square deviation (RMSD) of the natural ligand binding posture was used as the original RMSD label data, with a value of 0.

[0032] The experimentally determined binding strength pkd was used as the original pkd tag data;

[0033] The binding strength PKD tag data of the virtual bound pose is calculated according to the following formula:

[0034] pkd label =pkd nat -W*RMSD real

[0035] In the above formula, pkd label To train a set of virtual combined pose PKD labeled data; PKD nat The binding strength pkd represents the natural binding attitude of the ligand; W is the decay factor of the binding strength pkd as RMSD increases; RMSD real The root mean square deviation of the virtual fused pose from the natural fused pose.

[0036] Furthermore, in step S2, the heterogeneous graph G RR,LL,RL Only protein atoms within the first distance threshold D1 of the ligand atom are considered.

[0037] Furthermore, the heterogeneous graph G RR,LL,RL In this context, the ligand nodes and the edges formed between them are defined as the ligand subgraph G. L ligand subgraph Protein nodes and the edges formed between them are defined as protein subgraphs G. R That is, protein subgraph

[0038] In ligand subgraph G LThe code defines seven ligand atom types: C, N, O, P, S, Hal, and DU; where Hal represents the halogen elements F, Cl, Br, and I; and DU represents element types other than C, N, O, P, S, and Hal; if a chemical bond exists between atoms i and j in the ligand, then the node... The edges are formed between The physicochemical information of the chemical bonds is extracted and used together with the side lengths to form the ligand subgraph G. L Edge features;

[0039] In protein subgraph G R In this study, protein atom nodes are characterized by the type of residue, element type, whether they belong to the main chain or side chain, aromaticity, charge, and distance to α-C. When two protein nodes... An edge is formed when the distance between them is less than the second distance threshold d2. side The length of the protein subgraph G R The edge features, the second distance threshold d2 = 5-10 angstroms;

[0040] If protein nodes With ligand nodes When the distance between them is less than the third distance threshold d3, a directed edge is formed. in, The direction is from arrive On the contrary, Indicates from arrive The third distance threshold d3 = 8-15 angstroms and d3 ≤ D1 ≤ 15 angstroms.

[0041] Furthermore, the heterogeneous graph G RR,LL,RL It also includes directional features used to describe the relative positions of protein-ligand atoms, specifically:

[0042] Introducing a geometric center, protein node, based on ligand binding posture. ligand nodes The dihedral angle formed by the four points of the α-C atom And introduce the included angles of two adjacent sides as θ1 and θ2 respectively;

[0043] The directional features include: cos(θ1), cos(θ2), and the distance from the protein atom to the ligand atom.

[0044] Furthermore, the heterogeneous graph G RR,LL,RL It also includes an undirected graph G describing the protein-binding pocket.p :

[0045]

[0046] In the above formula, This represents the nodes in the pocket diagram. Represents the edges formed between nodes;

[0047] The undirected graph G of the protein binding pocket p It contains the geometric features of the residues at the protein-binding pocket;

[0048] The geometric features include the following definitions:

[0049] Residues of the protein within a fourth distance threshold d4 of the ligand crystal structure constitute a binding pocket; wherein the fourth distance threshold d4 = 8-12 Å and d4 ≤ D1 ≤ 15 Å;

[0050] node This indicates residues located in the protein-binding pocket;

[0051] An edge is formed only when the distance between the α-C atoms i and j of two residues is less than the fifth distance threshold d5. Among them, the fifth distance threshold d5 = 8-12 angstroms;

[0052] Meanwhile, residue type, distance distribution between internal atoms, dihedral angle of main chain / side chain, and position relative to the pocket center are used as node features; the distance between main chain atoms between nodes is used as edge features.

[0053] Furthermore, step S4 specifically includes:

[0054] S401: Input the two graphs encoded in step S2 into the scoring model constructed in step S3, learn the graph representation through the feature extraction module in the edge-based graph attention network, and encode it into a protein-ligand interaction latent space of length 1024.

[0055] Step S402: The latent space is decoded by the RMSD decoder and pkd decoder to obtain the predicted root mean square deviation of ligand binding posture (RMSD) and protein-ligand binding strength (pkd), respectively. At the same time, the pkd decoder also predicts the decay factor W of the binding strength pkd as the root mean square deviation (RMSD) decreases.

[0056] Step S403: Substitute the predicted RMSD, pkd, and W from step S402 into the loss function, differentiate the loss function, and backpropagate the gradient layer by layer to update the neural network parameters in the direction of minimizing the loss function. This makes the RMSD and pkd predicted by the scoring model close to the RMSD and pkd label data in step S1, and finally obtains the trained scoring model.

[0057] Furthermore, in step S403, the loss function used for training the scoring model is:

[0058]

[0059] In the above formula, α, β, and γ are the summation weights; mse(rmsd pred ,rmsd real ) represents the mean squared error between the predicted RMSD and the true RMSD; mse(pkd) label ,pkd pred The mean squared difference between the predicted binding strength pkd and the pkd tag data; The correction term represents the average value of the predicted binding strength pkd in each batch during training. The purpose of this correction term is to control the convergence direction of the model.

[0060] Meanwhile, this invention also provides a method for optimizing the conformation of a protein pocket side chain-ligand. Based on the ligand binding posture predicted by the protein-ligand interaction prediction method described in any of the preceding claims, the conformation of the protein pocket side chain-ligand is optimized. The optimization method specifically includes:

[0061] Step S(1): First, encode the ligands into a 6+k dimensional vector; for N side chain residues in the protein pocket, encode each side chain i residue into x. i 3D vector, x i The length of is equal to the number of rotatable bonds in the side chain of residue i;

[0062] Step S(2): Through the 6+k dimensional vector of the ligand and the N x-values ​​of the protein pocket side chain residues i 3D vectors are used to reconstruct the three-dimensional structure of ligands and proteins;

[0063] Step S(3): Using the scoring model trained in the protein-ligand interaction prediction method, predict the root mean square deviation of ligand binding posture (RMSD) and the protein-ligand binding strength (pkd).

[0064] The interatomic distances within the ligand, the interatomic distances between the protein and ligand, and the interatomic distances between the side chains of the protein pocket are calculated and used as geometric constraint terms.

[0065] Integrate RMSD, pkd, and geometric constraint terms to obtain the final comprehensive score for RMSD, pkd, and geometric constraint terms;

[0066] Step S(4): Take the gradient derivative of the comprehensive score in step S(3), and return the gradient to the 6+k dimensional vector of the ligand and the x-value of each residue of the pocket side chain. i The dimension vector is used to obtain the latest representation vector;

[0067] Step S(5): Iterate through steps S(2)-S(4) until the overall score no longer decreases, and then obtain the final optimized pocket side chain-ligand conformation.

[0068] Further, in step S(3), the geometric constraint term includes:

[0069] The formula for calculating the interatomic distance within a ligand is:

[0070]

[0071] The formula for calculating the distance between protein and ligand atoms is:

[0072]

[0073] The formula for calculating the interatomic distance between the side chains of a protein pocket is:

[0074]

[0075] In the above formula, σ lig σ pro-lig σ sidechain These are the interatomic distances within the ligand, the interatomic distances between the protein and ligand, and the interatomic distances between the side chains of the protein pocket; N lig N pro N ha These are the number of ligand atoms, the number of protein atoms interacting with the ligand, and the number of protein pocket side chain atoms, respectively. These are the cutoff values ​​of e. lig -d ij cutpff pro-lig -d ij cutoff sidechain -d ij Power; d ij The distance between two atoms i and j; cutpff pro-lig cutoff sidechainThese are the distance thresholds for the interaction between two atoms, namely, atoms inside the ligand, atoms between the protein and the ligand, and atoms between the side chains of the protein pocket. Repulsion only occurs when the distance between the two atoms is less than the corresponding distance threshold.

[0076] Compared with the prior art, the beneficial effects of the present invention are:

[0077] (1) The scoring model of the present invention uses two graphs to characterize protein-ligand interactions, namely the protein activity pocket graph and the protein-ligand atom interaction graph. The latent space of protein-ligand interactions is encoded by an edge-based graph attention network. The latent space is decoded by two decoders to obtain the root mean square deviation of ligand binding posture (RMSD) and the binding strength (pkd), respectively. This enables the scoring model to simultaneously predict the root mean square deviation of ligand binding posture (RMSD) and the binding strength (pkd). At the same time, the scoring model of the present invention can have excellent performance in scoring, ranking, docking and screening tasks, which improves the universality and transferability of models in this field.

[0078] (2) In the training process of the scoring model of the present invention, a decay factor W is introduced to predict the binding strength as RMSD increases, and based on this, the binding strength pkd label data of virtual binding pose is defined, and a new loss function is constructed. By taking the derivative of the loss function, backpropagating and updating the model parameters, the final trained scoring model can achieve good prediction accuracy and can achieve balanced performance in scoring, sorting, docking and screening tasks, overcoming the technical problem that existing scoring models cannot achieve satisfactory results in these four tasks at the same time.

[0079] (3) The protein-ligand optimization method proposed in this invention can simultaneously optimize the protein pocket side chain and ligand structure, thereby more accurately restoring the specific binding posture of the ligand in the protein pocket. That is, the deep learning model is realized in the dynamics of protein-ligand, overcoming the technical difficulty of realizing the fine modeling of protein-ligand conformation in the current fine modeling of protein binding pocket.

[0080] (4) Since the scoring model of the present invention has high universality and transferability in the field, in practical applications, the scoring model and optimization method of the present invention can realize the functions of various models, such as: ligand binding posture assessment, protein-ligand binding strength assessment, and fine modeling of protein-ligand complexes. The binding strength value given by the model corresponds to the binding energy measured in experiments. Therefore, it is still applicable to the prediction of two-body interactions in other scenarios, such as protein-protein interactions and protein-peptide interactions. Its application value and theoretical guidance significance are significant. Attached Figure Description

[0081] Figure 1 This is a flowchart of the scoring model algorithm in an embodiment of the present invention;

[0082] Figure 2 This is a schematic diagram illustrating the directional features of the relative positions of protein-ligand atoms in an embodiment of the present invention;

[0083] Figure 3 This is a flowchart of the protein-ligand optimization framework algorithm in an embodiment of the present invention;

[0084] Figure 4 Figure 1 shows a comparison of the ligand structure before and after optimization in an embodiment of the present invention. Figure 2 shows the ligand structure before optimization and Figure 3 shows the ligand structure after optimization. Detailed Implementation

[0085] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0086] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this disclosure belongs. The terminology used herein is for the purpose of describing embodiments of the invention only and is not intended to limit the invention.

[0087] For ease of description, the following explanations are provided regarding the terminology that may be involved in the method of this invention:

[0088] Docking Pose: The ligand conformation, or binding posture, produced by molecular docking.

[0089] RMSD: Root Mean Square Deviation, used to characterize the conformational difference between the ligand Docking Pose and its native state.

[0090] pkd: Binding strength, representing how tightly the ligand binds in the protein binding pocket.

[0091] Scoring, Ranking, Docking, and Screening Power: These refer to the scoring, ranking, matching, and screening capabilities, respectively. These are four metrics for evaluating a scoring model, each corresponding to different application scenarios.

[0092] Example 1

[0093] Combination Figure 1As shown in the figure, this invention provides a method for predicting protein-ligand interactions, the method comprising the following steps:

[0094] Step S1: Construct the dataset and create labeled data;

[0095] Using the 2019 version of the PDBbind database (PDBbind v2019) and experimentally determined binding strengths as the original dataset, tagged data was created. This tagged data included RMSD tag data and pkd tag data.

[0096] Specifically, the root mean square deviation (RMSD) of the binding posture of these natural ligands is used as RMSD tag data, with a value of 0; the binding strength (pkd) determined experimentally is used as pkd tag data.

[0097] Wherein, the binding strength pkd is the negative logarithm of the dissociation constant KD / inhibition constant KI / half-inhibition concentration IC50, that is:

[0098] pkd = -log(KD) or -log(KI) or -log(IC50)

[0099] Meanwhile, PDBbind v2019 includes a subset called the refine set, from which 1000 natural protein-ligand complex samples are randomly selected to form the validation set;

[0100] PDBbind v2019 includes a subset called the core set (CASF-2016core set), which contains 285 protein-ligand complex samples that make up the test set;

[0101] The training set in PDBbind v2019 consists of the samples remaining after excluding the validation set (1000 samples) and the test set (285 samples).

[0102] Furthermore, molecular docking can be used to generate more virtual binding poses for ligands in the training and validation sets, thereby further expanding the sample data. The RMSD of these virtual binding poses relative to the natural binding poses can be directly calculated, but the binding strength between the virtual binding poses and proteins cannot be directly calculated. Therefore, this invention proposes an assumption: the natural binding pose of the ligand in the protein has the greatest binding strength, while the binding strength pkd of the virtual binding poses (poses) generated by molecular docking decreases with increasing RMSD.

[0103] Therefore, based on the above assumptions, a decay factor W is introduced to reduce the binding strength pkd as the RMSD increases. Based on the binding strength pkd of the ligand’s natural binding posture, the calculated RMSD, and the decay factor W, the pkd tag data of the ligand’s virtual binding posture is calculated.

[0104] The following formula is given:

[0105] pkd label =pkd nat -W*RMSD real (1)

[0106] In equation (1) above, pkd label PKD tag data for ligand virtual binding posture; PKD nat ρkd represents the binding strength of the ligand in its natural binding posture; W is the decay factor of the binding strength pkd as RMSD increases; RMSD real This represents the root mean square deviation of the ligand binding posture relative to the natural binding posture.

[0107] At this point, the label data (RMSD label data and pkd label data) of the samples in both the training and validation sets can be obtained.

[0108] Step S2: Establish a protein-ligand complex model;

[0109] The protein-ligand relationship is encoded into two graphs: a protein binding pocket graph and a protein-ligand atom interaction graph.

[0110] The protein-ligand atomic interaction diagram is an isomer diagram of protein-ligand interactions, G. RR,LL,RL The heterogeneous graph G RR,LL,RL The model does not model all the atoms of the protein, but only considers the protein atoms within the first distance threshold D1 of the ligand atom. The first distance threshold D1 can be set according to the actual situation, generally not exceeding 15 angstroms, preferably 8-12 angstroms; in the embodiment of the present invention, it is specifically taken as 8 angstroms.

[0111] The heterogeneous graph G RR,LL,RL There are two types of nodes and four types of edges in the system.

[0112] The two types of nodes are: protein atom nodes. and ligand atomic nodes The diagram uses circles to represent these areas.

[0113] The four types of edges are: white matter node - protein node. Ligand node - ligand node Protein node-ligand node Ligand nodes - protein nodes Meanwhile, these four types of edges are represented by two line types in the graph: dashed lines represent edges between protein nodes and between ligand nodes, while straight lines represent edges between protein nodes and between ligand nodes and between protein nodes.

[0114] Heterogeneous graph G RR,LL,RL In this context, the ligand nodes and the edges formed between them are defined as the ligand subgraph G. L ligand subgraph Meanwhile, in the ligand subgraph G L The standard defines seven ligand atom types: C, N, O, P, S, Hal, and DU, where Hal represents the halogen elements F, Cl, Br, and I, and DU represents element types other than C, N, O, P, S, and Hal. If a chemical bond exists between atoms i and j in the ligand, then the node... The edges are formed between Furthermore, the physicochemical information of chemical bonds was extracted using the rdkit package in Python, and this information, along with the edge length, was used as an edge feature.

[0115] Heterogeneous graph G RR,LL,RL In this context, protein nodes and the edges formed between them are defined as protein subgraphs G. R That is, protein subgraph Meanwhile, in protein subgraph G R In this context, protein atom nodes are characterized by factors such as the type of residue they occupy, element type, whether they belong to the main chain or side chain, aromaticity, charge, and distance to the α-C atom (the carbon atom directly connected to the functional group). When two protein nodes... An edge is formed when the distance is less than the second distance threshold d2. side The length of the edge is used as the edge feature, where the second distance threshold d2 = 5 angstroms.

[0116] If protein nodes With ligand nodes When the distance between them is less than the third distance threshold d3, a directed edge is formed. in The direction is from arrive On the contrary, Indicates from arrive The third distance threshold d3 = 8 angstroms.

[0117] It should be noted that in some other embodiments, the second distance threshold d2 and the third distance threshold d3 can also be adjusted according to the actual situation. Specifically, the adjustment range of the second distance threshold d2 is not limited by the first distance threshold D1, and is preferably 5-10 angstroms; the adjustment range of the third distance threshold d3 must be less than or equal to the first distance threshold D1, and is 8-15 angstroms, preferably 8-12 angstroms.

[0118] To more comprehensively describe the relative positions between protein-ligand atoms, in addition to applying interatomic distances, embodiments of the present invention also propose directional features for describing the relative positions of protein-ligand atoms.

[0119] Specifically, a geometric center and protein node based on the ligand binding posture are introduced. ligand nodes The dihedral angle formed by the four points of the α-C atom And introduce the included angles between two adjacent sides as θ1 and θ2; for example Figure 2 As shown in the figure, P is the geometric center of the ligand binding posture, L is an atom of the ligand binding posture, R represents an atom on the residue, Cα represents the α-C of the protein residue, and the dihedral angle is... Let θ1 be the dihedral angle formed by points L, R, Cα, and α-C, and let θ2 be the angle between side Cα-R and RL. Finally, cos(θ1), cos(θ2), and the distance from the protein atom to the ligand atom are used together as features to describe the protein-ligand atom interaction edge (i.e., the introduced directional features).

[0120] Although the aforementioned isomorphic graph contains protein atomic information, the atomic-level node representation cannot directly describe the internal state of residues. Therefore, this embodiment of the invention also constructs an undirected graph describing the protein binding pocket. in, This represents the nodes in the pocket diagram. This represents the edges formed between nodes.

[0121] The undirected graph G of the protein binding pocket p This includes the geometric features of residues at the protein-binding pocket, defined as follows: residues of the protein located less than a fourth distance threshold d4 from the ligand crystal structure constitute the binding pocket, where the fourth distance threshold d4 = 8 Å. This refers to residues at the protein binding pocket. An edge is formed only when the distance between the α-C atoms i and j of two residues is less than the fifth distance threshold d5. The fifth distance threshold d5 = 10 angstroms; at the same time, the residue type (represented by one-heat encoding), the distance distribution between internal atoms, the dihedral angle of the main / side chain, and the position relative to the pocket center are used as node features, and the distance between main chain atoms between nodes is used as edge features.

[0122] It should be noted that in some other embodiments, the fourth distance threshold d4 and the fifth distance threshold d5 can also be adjusted according to the actual situation. Specifically, the adjusted value of the fourth distance threshold d4 must be less than or equal to the first distance threshold D1, and the adjustment range is 8-12 angstroms; the adjusted value of the fifth distance threshold d5 is not limited by the first distance threshold D1, and the adjustment range is also 8-12 angstroms.

[0123] The directional features of the relative positions of protein-ligand atoms described above can be understood as a refinement of the protein-ligand complex model. Because nodes and edges in a graph have their own characteristics, the directional information in the protein-ligand graph is actually one of the features of protein-ligand edges and ligand-protein edges. The undirected graph of protein binding pockets can be understood as a supplement to the interaction description, because the binding of proteins and ligands is affected by the environment and physicochemical properties of the binding pocket. The purpose of constructing an undirected graph of protein binding pockets is to integrate the information of the protein binding pockets into the characterization of protein-ligand interactions.

[0124] Step S3: Construct an edge-based graph attention network (EdgeGAT layer) scoring model;

[0125] The scoring model architecture consists of three parts: a feature extraction module, an RMSD decoder, and a pkd decoder.

[0126] The feature extraction module includes two independent edge-based graph attention layers for encoding features from the protein binding pocket graph and the protein-ligand atom interaction graph into a latent space of protein-ligand interactions with a length of 1024.

[0127] The RMSD decoder and the pkd decoder are used to decode the latent space to obtain the predicted root mean square deviation of ligand binding attitude (RMSD) and protein-ligand binding strength (pkd), respectively.

[0128] Step S4: Train the scoring model;

[0129] The two graphs encoded in step S2 are input into the scoring model (EdgeGAT layer based graph attention network) constructed in step S3. After feature extraction, encoding and decoding of the protein-ligand interaction latent space, backpropagation and parameter update, the trained scoring model is obtained.

[0130] Step S4 specifically includes:

[0131] S401: Input the two graphs encoded in step S2 into the scoring model (EdgeGAT layer based graph attention network) constructed in step S3, learn the graph representation through the feature extraction module in the edge-based graph attention network, and encode it into a protein-ligand interaction latent space of length 1024.

[0132] Step S402: The latent space is decoded by the RMSD decoder and pkd decoder to obtain the predicted root mean square deviation of ligand binding posture (RMSD) and protein-ligand binding strength (pkd), respectively. At the same time, the pkd decoder also predicts the decay factor W of the binding strength pkd as the root mean square deviation (RMSD) decreases.

[0133] Step S403: Substitute the predicted RMSD, pkd, and W from step S402 into the loss function, differentiate the loss function, and backpropagate the gradient layer by layer to update the neural network parameters in the direction of minimizing the loss function. This makes the RMSD and pkd predicted by the scoring model close to the RMSD and pkd label data in step S1, and finally obtains the trained scoring model.

[0134] During the training process, the binding strength pkd label data of the virtual combined pose in the training set is also calculated according to formula (1) in step S1.

[0135] The loss function used for training the scoring model is:

[0136]

[0137] In equation (2) above, α, β, and γ are the summation weights; mse(rmsd pred ,rmsd real ) represents the mean squared error between the predicted RMSD and the true RMSD; mse(pkd) label ,pkd pred The mean square error between the predicted binding strength pkd and the pkd tag is denoted as . The correction term represents the average value of the predicted binding strength pkd in each batch during training. The purpose of this correction term is to control the convergence direction of the model.

[0138] Step S5: Using the trained scoring model, predict the interaction between the protein and ligand to be tested, and obtain the predicted root mean square deviation of ligand binding posture (RMSD) and protein-ligand binding strength (pkd).

[0139] Example 2

[0140] This invention provides a method for optimizing the conformation of protein pocket side chain-ligand, which optimizes the conformation of protein pocket side chain-ligand based on the ligand binding posture predicted by the protein-ligand interaction prediction method in Example 1.

[0141] like Figure 3 As shown, the optimization method specifically includes:

[0142] Step S(1): First, encode the ligands into a 6+k dimensional vector; for N side chain residues in the protein pocket, encode each side chain i residue into x. i 3D vector, x i The length of the vector is equal to the number of rotatable bonds in the side chain of residue i;

[0143] Step S(2): Through the 6+k dimensional vector of the ligand and the N x-values ​​of the protein pocket side chain residues i 3D vectors are used to reconstruct the three-dimensional structure of ligands and proteins;

[0144] Step S(3) uses the scoring model trained in Example 1 to predict the root mean square deviation of ligand binding posture (RMSD) and protein-ligand binding strength (pkd).

[0145] The interatomic distances within the ligand, the protein-ligand interatomic distances, and the interatomic distances between protein pocket side chains are calculated and used as geometric constraint terms; where:

[0146] The formula for calculating the interatomic distance within a ligand is:

[0147]

[0148] The formula for calculating the distance between protein and ligand atoms is:

[0149]

[0150] The formula for calculating the interatomic distance between the side chains of a protein pocket is:

[0151]

[0152] In equations (3)-(5) above, σ lig σ pro-lig σ sidechain These are the interatomic distances within the ligand, the interatomic distances between the protein and ligand, and the interatomic distances between the side chains of the protein pocket; N lig N pro N ha These are the number of ligand atoms, the number of protein atoms interacting with the ligand, and the number of protein pocket side chain atoms, respectively. These are the cutoff values ​​of e. lig -d ij cutpff pro-lig -d ij cutoff sidechain -d ij Power; d ij The distance between two atoms i and j; cutpff pro-lig cutoff sidechain These represent the distance thresholds for two-atom interactions, namely, interactions between atoms within the ligand, between protein and ligand atoms, and between atoms in the protein pocket side chain. Repulsive interactions only occur when the distance between the two atoms is less than the corresponding distance threshold. cutpff pro-lig cutoff sidechain Typically, the value is between 2.5 and 3.5 angstroms, and the specific value of each distance threshold can be set according to the actual situation.

[0153] Integrate RMSD, pkd, and geometric constraint terms to obtain the final comprehensive score for RMSD, pkd, and geometric constraint terms;

[0154] Step S(4): Take the gradient derivative of the comprehensive score in step S(3), and return the gradient to the 6+k dimensional vector of the ligand and the x-value of each residue of the pocket side chain. i The dimension vector is used to obtain the latest representation vector;

[0155] Step S(5): Iterate through steps S(2)-S(4) until the overall score no longer decreases, and then obtain the final optimized protein pocket side chain-ligand conformation.

[0156] Example 3

[0157] The present invention verifies the methods in Embodiment 1 and Embodiment 2 respectively, and the verification results are as follows:

[0158] 1. Gain effect of scoring model:

[0159] The scoring model in Example 1 was evaluated on the CASF-2016 standard, a mainstream evaluation standard for scoring functions, and compared with the latest mainstream scoring functions. The results are shown in Tables 2-4 below:

[0160] Table 2. Performance of different scoring models in the Scoring and Ranking tasks

[0161] Model Name Scoring ability (PCC) Sorting Capability (SCC) pkd-model 0.825 0.707 DeepRMSD+Vina 0.547 0.49 RTMScore 0.455 0.529

[0162] The evaluation metrics for the Scoring and Ranking tasks are the Pearson Correlation Coefficient (PCC) and the Spearman Correlation Coefficient (SCC), respectively.

[0163] As can be seen from Table 2, in the Scoring and Ranking tasks, the scoring model (RMSD-model and pkd-model) in Embodiment 1 of this invention is significantly better than the other two models.

[0164] Table 3. Performance of different scoring models in the Docking task

[0165]

[0166] As can be seen from Table 3, in the docking task, the RMSD prediction module (RMSD-model) of the scoring model in Embodiment 1 of the present invention outperforms other current scoring models in both cases where natural ligand binding posture is included and where natural ligand binding posture is not included. Moreover, the pkd prediction module (pkd-model) of the scoring model in Embodiment 1 of the present invention still outperforms other scoring models in the Top 1 success rate when natural binding posture is not included.

[0167] Table 4. Performance of different scoring models in the Screening Power task

[0168] Model Name Enrichment factor (1%) Top 1 success rate pkd-model 18.36 70.2％ RTMScore 28.0 66.7％ DeepRMSD+Vina 21.95 47.4％

[0169] In the Screening Power task, two evaluation metrics were used to measure the screening ability of the scoring model: enrichment factor (1%) and Top 1 success rate. As shown in Table 4, the pkd-model of Example 1 of this invention performs exceptionally well in identifying ligands with the highest affinity, outperforming other models.

[0170] Furthermore, the protein-ligand conformation optimization method described in Example 2 of this invention was used to optimize the protein-ligand conformation, and the results are as follows. Figure 4 As shown in the figure, red (color not shown) represents the native conformation, and blue (color not shown) represents the virtual conformation of molecular docking.

[0171] from Figure 4As can be seen, the pseudoconformations generated by molecular docking have the ability to be reduced to a near-native state, demonstrating the reliability and practicality of the protein pocket side chain-ligand conformation optimization method in Example 2 of this invention.

[0172] The above description is merely an embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the scope of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for predicting protein-ligand interactions, characterized in that, The method includes the following steps: Step S1: Construct the dataset and create labeled data; Step S2: Establish a protein-ligand complex model; The protein-ligand relationship is encoded into two diagrams: a protein binding pocket diagram and a protein-ligand atom interaction diagram. The protein-ligand atomic interaction diagram is an isomer diagram of protein-ligand interactions, G. RR,LL,RL The heterogeneous graph G RR,LL,RL There are two types of nodes and four types of edges in the system; The two types of nodes are: protein atom nodes. and ligand atomic nodes The four types of edges are: protein node - protein node Ligand node - ligand node Protein node-ligand node Ligand nodes - protein nodes Step S3: Construct an edge-based graph attention network scoring model; The scoring model includes: a feature extraction module, an RMSD decoder, and a pkd decoder; The feature extraction module includes two independent edge-based graph attention layers for encoding features from the protein binding pocket graph and the protein-ligand atom interaction graph into a latent space for protein-ligand interactions. The RMSD decoder and the pkd decoder are used to decode the latent space to obtain the predicted root mean square deviation of ligand binding attitude (RMSD) and protein-ligand binding strength (pkd), respectively. Step S4: Train the scoring model; The two graphs encoded in step S2 are input into the scoring model constructed in step S3. After feature extraction, encoding and decoding of the protein-ligand interaction latent space, backpropagation and parameter update, the trained scoring model is obtained. Step S5: Using the trained scoring model, predict the interaction between the protein and ligand to be tested, and obtain the predicted root mean square deviation of ligand binding posture (RMSD) and protein-ligand binding strength (pkd).

2. The prediction method according to claim 1, characterized in that, The dataset includes a validation set, a test set, and a training set. More virtual binding postures are generated for the ligands in the training and validation sets through molecular docking, thereby expanding the sample data of the training and validation sets. The dataset includes RMSD label data and pkd label data; The root mean square deviation (RMSD) of the natural ligand binding posture was used as the original RMSD label data, with a value of 0. The experimentally determined binding strength pkd was used as the original pkd tag data; The binding strength PKD tag data of the virtual bound pose is calculated according to the following formula: pkd label =pkd nat -W*RMSD real In the above formula, pkd label To train a set of virtual combined pose PKD labeled data; PKD nat The binding strength pkd represents the natural binding attitude of the ligand; W is the decay factor of the binding strength pkd as RMSD increases; RMSD real The root mean square deviation of the virtual fused pose from the natural fused pose.

3. The prediction method according to claim 1, characterized in that, In step S2, the heterogeneous graph G RR,LL,RL Only protein atoms within the first distance threshold D1 of the ligand atom are considered.

4. The prediction method according to claim 3, characterized in that, The heterogeneous graph G RR,LL,RL In this context, the ligand nodes and the edges formed between them are defined as the ligand subgraph G. L ligand subgraph Protein nodes and the edges formed between them are defined as protein subgraphs G. R That is, protein subgraph In ligand subgraph G L The code defines seven ligand atom types: C, N, O, P, S, Hal, and DU; where Hal represents the halogen elements F, Cl, Br, and I; and DU represents element types other than C, N, O, P, S, and Hal; if a chemical bond exists between atoms i and j in the ligand, then the node... The edges are formed between The physicochemical information of the chemical bonds is extracted and used together with the side lengths to form the ligand subgraph G. L Edge features; In protein subgraph G R In this study, protein atom nodes are characterized by the type of residue, element type, whether they belong to the main chain or side chain, aromaticity, charge, and distance to the α-C atom. When two protein nodes... An edge is formed when the distance between them is less than the second distance threshold d2. side The length of the protein subgraph G R The edge features, the second distance threshold d2 = 5-10 angstroms; If protein nodes With ligand nodes When the distance between them is less than the third distance threshold d3, a directed edge is formed. in, The direction is from arrive On the contrary, Indicates from arrive The third distance threshold d3 = 8-15 angstroms and d3 ≤ D1 ≤ 15 angstroms.

5. The prediction method according to claim 3, characterized in that, The heterogeneous graph G RR,LL,RL It also includes directional features used to describe the relative positions of protein-ligand atoms, specifically: Introducing a geometric center, protein node, based on ligand binding posture. ligand nodes The dihedral angle formed by the four points of the α-C atom And introduce the included angles of two adjacent sides as θ1 and θ2 respectively; The directional features include: cos(θ1), cos(θ2), and the distance from the protein atom to the ligand atom.

6. The prediction method according to claim 3, characterized in that, The heterogeneous graph G RR,LL,RL It also includes an undirected graph G describing the protein-binding pocket. p : In the above formula, ε represents the node in the combined pocket graph. p Represents the edges formed between nodes; The undirected graph G of the protein binding pocket p It contains the geometric features of the residues at the protein-binding pocket; The geometric features include the following definitions: Residues of the protein within a fourth distance threshold d4 of the ligand crystal structure constitute a binding pocket; wherein the fourth distance threshold d4 = 8-12 Å and d4 ≤ D1 ≤ 15 Å; node This indicates residues located in the protein-binding pocket; An edge is formed only when the distance between the α-C atoms i and j of two residues is less than the fifth distance threshold d5. Among them, the fifth distance threshold d5 = 8-12 angstroms; Meanwhile, residue type, distance distribution between internal atoms, dihedral angle of main chain / side chain, and position relative to the pocket center are used as node features; the distance between main chain atoms between nodes is used as edge features.

7. The prediction method according to claim 1, characterized in that, Step S4 specifically includes: S401: Input the two graphs encoded in step S2 into the scoring model constructed in step S3, learn the graph representation through the feature extraction module in the edge-based graph attention network, and encode it into a protein-ligand interaction latent space of length 1024. Step S402: The latent space is decoded by the RMSD decoder and pkd decoder to obtain the predicted root mean square deviation of ligand binding posture (RMSD) and protein-ligand binding strength (pkd), respectively. At the same time, the pkd decoder also predicts the decay factor W of the binding strength pkd as the root mean square deviation (RMSD) decreases. Step S403: Substitute the predicted RMSD, pkd, and W from step S402 into the loss function, differentiate the loss function, and backpropagate the gradient layer to update the neural network parameters in the direction of minimizing the loss function. This makes the RMSD and pkd predicted by the scoring model close to the RMSD and pkd label data in step S1, and finally obtains the trained scoring model.

8. The prediction method according to claim 7, characterized in that, In step S403, the loss function used for training the scoring model is: In the above formula, α, β, and γ are the summation weights; mse(rmsd pred ,rmsd real ) represents the mean squared error between the predicted RMSD and the true RMSD; mse(pkd) label ,pkd pred The mean squared difference between the predicted binding strength pkd and the pkd tag data; The correction term represents the average value of the predicted binding strength pkd in each batch during training. The purpose of this correction term is to control the convergence direction of the model.

9. A method for optimizing the conformation of a protein pocket side chain-ligand, comprising optimizing the conformation of the protein pocket side chain-ligand based on the ligand binding posture predicted by the protein-ligand interaction prediction method as described in any one of claims 1-8, characterized in that, The optimization method specifically includes: Step S(1): First, encode the ligands into a 6+k dimensional vector; for N side chain residues in the protein pocket, encode each side chain i residue into x. i 3D vector, x i The length of is equal to the number of rotatable bonds in the side chain of residue i; Step S(2): Through the 6+k dimensional vector of the ligand and the N x-values ​​of the protein pocket side chain residues i 3D vectors are used to reconstruct the three-dimensional structure of ligands and proteins; Step S(3): Using the scoring model trained in the protein-ligand interaction prediction method, predict the root mean square deviation of ligand binding posture (RMSD) and the protein-ligand binding strength (pkd). The interatomic distances within the ligand, the interatomic distances between the protein and ligand, and the interatomic distances between the side chains of the protein pocket are calculated and used as geometric constraint terms. Integrate RMSD, pkd, and geometric constraint terms to obtain the final comprehensive score for RMSD, pkd, and geometric constraint terms; Step S(4): Take the gradient derivative of the comprehensive score in step S(3), and return the gradient to the 6+k dimensional vector of the ligand and the x-value of each residue of the pocket side chain. i The dimension vector is used to obtain the latest representation vector; Step S(5): Iterate through steps S(2)-S(4) until the overall score no longer decreases, and then obtain the final optimized pocket side chain-ligand conformation.

10. The optimization method according to claim 9, characterized in that, In step S(3), the geometric constraint term includes: The formula for calculating the interatomic distance within a ligand is: The formula for calculating the distance between protein and ligand atoms is: The formula for calculating the interatomic distance between the side chains of a protein pocket is: In the above formula, σ lig σ pro-lig σ sidechain These are the interatomic distances within the ligand, the interatomic distances between the protein and ligand, and the interatomic distances between the side chains of the protein pocket; N lig N pro N ha These are the number of ligand atoms, the number of protein atoms interacting with the ligand, and the number of protein pocket side chain atoms, respectively. These are the cutoff values ​​of e. lig -d ij cutpff pro-lig -d ij cutoff sidechain -d ij Power; d ij The distance between two atoms i and j; cutpff pro-lig cutoff sidechain These are the distance thresholds for the interaction between two atoms, namely, atoms inside the ligand, atoms between the protein and the ligand, and atoms between the side chains of the protein pocket. Repulsion only occurs when the distance between the two atoms is less than the corresponding distance threshold.

Images