A drug molecule virtual screening method based on binding mode analysis

Abstract

The present invention discloses a virtual screening method for drug molecules based on binding mode analysis, which belongs to the field of artificial intelligence drug discovery. PLIP is used to calculate the binding mode between the target to be screened and the small molecule, and the binding mode is encoded. The DBSCAN algorithm is used to cluster the binding modes to identify groups of compounds with similar binding modes. In the screening stage, two key indicators are introduced: an energy index is formed based on the hydrogen bond compensation number and the AutoDock Vina score, and an activity index is used based on the correlation between the relevant residues in the binding mode and the alpha carbon atom of the key site, which is used to evaluate the binding stability and characteristic activity between the target and the small molecule, and to perform multi-target frontier analysis on the candidate molecules. The method of the present invention significantly improves the efficiency and accuracy of virtual screening, enhances the detailed analysis of intermolecular interactions, is suitable for small-scale and large-scale molecular screening, reduces the cost and time of drug research and development, enriches the diversity of virtual screening results, and provides more potential candidate molecules for subsequent experiments.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of artificial intelligence drug discovery, and particularly relates to a drug molecule virtual screening method based on binding mode analysis. BACKGROUND

[0002] Virtual screening technology has become an indispensable part of modern drug discovery since its inception. Through computer simulation, scientists can screen a large number of small molecule compounds in a short time, predict their binding affinity with target proteins, and thus provide candidate compounds for subsequent experimental verification. The core of this method is the computational analysis of the interaction between molecular structure and protein target, especially in the early stage of drug development, virtual screening can greatly reduce the screening time, reduce the cost, and improve the screening efficiency.

[0003] Traditional virtual screening methods usually rely on molecular docking programs such as AutoDockVina, which predict the binding conformation and binding energy of small molecules with target proteins to judge the activity of compounds. However, this method is mainly based on the level of docking score, which may ignore some key details of molecular interactions such as hydrogen bonds, hydrophobic interactions, salt bridges, etc. These details often determine the stability and specificity of the molecule binding to the target, and are the key to screening high-quality candidate molecules. Therefore, with the development of virtual screening technology, more and more researchers have begun to focus on binding mode analysis in order to more accurately identify active molecules and optimize screening results.

[0004] The introduction of binding mode analysis has brought a qualitative leap to virtual screening. By deeply analyzing the binding mode of small molecules and target proteins, scientists can not only identify key sites that contribute to binding stability, but also better understand the similarities and differences between different compounds. This analysis method not only improves the accuracy of screening, but also provides new ideas for drug design, such as enhancing certain key interactions through chemical modification to improve drug activity and selectivity. In addition, binding mode analysis can also help identify and classify compounds with similar binding modes, thereby enriching the diversity of virtual screening results and providing a wider candidate set for subsequent experimental screening.

[0005] Although binding mode analysis provides more depth and accuracy for virtual screening, relying solely on this technology is still insufficient to fully optimize the diversity and accuracy of screening results. For example, existing methods may not be able to fully assess the comprehensive impact of intermolecular interactions on binding stability and binding mode characteristic activity in some cases. Therefore, further technical improvements are needed to make up for these deficiencies in order to achieve more comprehensive and accurate screening results. SUMMARY

[0006] To solve the above problems, the application discloses a drug molecule virtual screening method based on binding mode analysis. The binding mode between the target to be screened and small molecules is calculated by using PLIP, the binding mode is encoded, and the DBSCAN algorithm is used for clustering the binding mode to identify the compound group with similar binding mode. In the screening stage, two key indicators are further introduced: the energy indicator is formed according to the hydrogen bond compensation number and the AutoDockVina score, and the activity indicator is formed according to the correlation between the relevant residues in the binding mode and the alpha carbon atoms of the key site, which is used to evaluate the binding stability between the target and the small molecule and the characteristic activity of the binding mode, and the multi-objective frontier analysis is performed on the candidate molecules.

[0007] To achieve the above purpose, the technical scheme of the application is as follows:

[0008] The application provides a drug molecule virtual screening method based on binding mode analysis, comprising the following steps:

[0009] Step S100: preparing a target to be screened and a small molecule data set which has been processed by AutoDockVina docking;

[0010] Step S200: using PLIP to calculate the binding mode between the target to be screened and small molecules, encoding the binding mode, and using the DBSCAN algorithm to cluster the distance matrix of the binding mode;

[0011] Step S300: selecting TopK1 molecules in each class, and calculating the hydrogen bond compensation number between the protein and the ligand;

[0012] Step S400: calculating the correlation between the residues involved in the binding mode of TopK1 molecules and the alpha carbon atoms of the key site;

[0013] Step S500: forming the energy indicator according to the hydrogen bond compensation number and the AutoDockVina score, and forming the activity indicator according to the correlation between the relevant residues in the binding mode and the alpha carbon atoms of the key site, performing multi-objective frontier analysis on the selected TopK1 molecules, and selecting the Pareto optimal solution set.

[0014] Further, step S100 specifically comprises the following steps:

[0015] A target protein to be screened is selected, and a relevant small molecule compound data set is prepared; the target and the small molecule data set are docked by using AutoDockVina, and the conformation file after docking is generated.

[0016] Further, step S200 specifically comprises the following steps:

[0017] Step S201: Merge the target protein and small molecule compound in step S100 into a complex using PyMOL, add hydrogen atoms to the small molecule, and delete excess water molecules in the complex.

[0018] Step S202: Calculate the binding mode of the complex in step S201 using PLIP, which includes 5 types of non-covalent bond interactions: Hydrogen Bonds, Hydrophobic Interactions, Salt Bridges, Pi-Cation Interactions, and Metal Complexes, as well as corresponding residue types, residue numbers, and chain names.

[0019] Step S203: Encode the binding mode between the target protein and small molecule complex calculated in step S202, and use the Tanimoto coefficient to measure similarity and calculate the distance matrix between binding modes.

[0020] Step S204: Cluster the distance matrix of binding modes in step S203 using the DBSCAN algorithm to identify groups of molecules with similar binding characteristics.

[0021] Among them, the above non-covalent bond interactions are explained as follows:

[0022] Hydrogen Bonds (Hydrogen Bonds): Hydrogen bonds are weak attractions between hydrogen atoms and highly electronegative atoms due to charge differences.

[0023] Hydrophobic Interactions (Hydrophobic Interactions): Hydrophobic interactions are non-covalent interactions that occur when hydrophobic molecules avoid water and aggregate in space.

[0024] Salt Bridges (Salt Bridges): Salt bridges are non-covalent electrostatic bonds between positively and negatively charged groups in proteins.

[0025] Pi-Cation Interactions (Pi-Cation Interactions): Pi-Cation Interactions are non-covalent attractions between π electrons on aromatic rings and cations.

[0026] Metal Complexes (Metal Complexes): Metal complexes are complexes formed by ligands and central metal ions through non-covalent coordination bonds, commonly found in interactions between proteins and metal-containing small molecules.

[0027] Further, step S300 specifically includes the following steps:

[0028] Step S301: according to the AutoDockVina score, top K1 molecules of each class in the clustering result in step S204 are selected, and a solvent-containing system file is generated by using the Gromacs command to process the target to be screened, the small molecule and the complex respectively.

[0029] Step S302: the number of hydrogen bonds in the solvent-containing system file generated in step S301 is calculated by using PyMOL and a set threshold, and the number of hydrogen bond compensations is calculated.

[0030] Further, step S400 specifically comprises the following steps:

[0031] Step S401: according to the AutoDockVina score, top K1 molecules of each class in the clustering result in step S204 are selected, and the target to be screened, the small molecule and the complex after docking processing are subjected to molecular dynamics simulation by using Gromacs, and a topology file and a trajectory file are generated;

[0032] Step S402: the topology file and the trajectory file generated in step S401 are loaded by using MDAnalysis, the Pearson correlation coefficient between the key site alpha carbon atom and the binding mode related residue obtained in step S202 is calculated, and then the correlation coefficients are added to obtain the total correlation coefficient.

[0033] Further, step S500 specifically comprises the following steps:

[0034] According to the hydrogen bond compensation number and the AutoDockVina score, an energy index is formed, and according to the correlation between the binding mode related residue and the key site alpha carbon atom, an activity index is formed, and multi-objective frontier analysis is performed on the selected top K1 molecules, and a Pareto optimal solution set is selected.

[0035] The beneficial effects of the present application are:

[0036] The application discloses a drug molecule virtual screening method based on binding mode analysis, and belongs to the field of artificial intelligence drug discovery. First, a target to be screened and a small molecule data set which have been processed by AutoDockVina docking are prepared; then, PLIP is used to calculate the binding mode between the target to be screened and the small molecules, the binding mode is encoded, and on this basis, the distance matrix of the binding mode is clustered by using a DBSCAN algorithm to identify a molecule group with similar binding characteristics; next, in the preliminary screening stage of the molecules, TopK1 molecules are selected from each class of clusters according to the AutoDockVina score, the number of hydrogen bond compensations between the target to be screened and the small molecules is further calculated, and the relevance of the residues involved in the binding mode of the TopK1 molecules and the alpha carbon atoms of the key sites is evaluated; finally, an energy index is formed according to the number of hydrogen bond compensations and the AutoDockVina score, and the relevance of the relevant residues in the binding mode and the alpha carbon atoms of the key sites is taken as an activity index, multi-objective frontier analysis is performed on the selected TopK1 molecules, and a Pareto optimal solution set is selected as the final screening result.

[0037] The method of the application systematically records the non-covalent interaction information in the binding mode through binding mode analysis, reduces redundant screening, and improves the efficiency and accuracy of virtual screening. At the same time, by comprehensively evaluating the number of hydrogen bond compensations and the relevance of the alpha carbon atoms of the key sites, the analysis of the interaction between molecules is further optimized. Ultimately, this systematic and refined analysis method, combined with multi-objective frontier analysis, not only significantly improves the effect of virtual screening, is suitable for small-scale and large-scale molecular screening, but also reduces the cost and time of drug research and development, enriches the diversity of virtual screening results, realizes more comprehensive and accurate screening results, and provides more potential candidate molecules for subsequent experiments. BRIEF DESCRIPTION OF DRAWINGS

[0038] Figure 1 The figure is a flowchart of the drug molecule virtual screening method based on binding mode analysis of the application.

[0039] Figure 2 The figure is a flowchart of the binding mode analysis calculation. DETAILED DESCRIPTION

[0040] The application will be further illustrated in combination with the drawings and specific embodiments. It should be understood that the following specific embodiments are only used to illustrate the application and are not used to limit the scope of the application.

[0041] The NMDA target is selected as an example target of the method.

[0042] The non-covalent bond interaction is explained as follows:

[0043] Hydrogen Bonds: Hydrogen bonds are weak attractive forces caused by the difference in charge between a hydrogen atom and a highly electronegative atom.

[0044] Hydrophobic Interactions: Hydrophobic interactions are non-covalent interactions that occur when hydrophobic molecules avoid water and aggregate in space.

[0045] Salt Bridges: Salt bridges are non-covalent electrostatic bonds formed between positively and negatively charged groups in proteins.

[0046] Pi-Cation Interactions(π-cation interactions): π-cation interactions are the non-covalent attractive forces between the π electrons on the aromatic ring and the cation.

[0047] Metal Complexes: Metal complexes are complexes formed by the binding of ligands to central metal ions through non-covalent coordination bonds, and are commonly seen in the interaction between proteins and metal-containing small molecules.

[0048] Example 1

[0049] Combine Figure 1 This embodiment discloses a virtual screening method for drug molecules based on binding mode analysis, which specifically includes the following steps:

[0050] Step S100: Prepare a dataset of the target NMDA to be screened and its 100W ligand small molecule that has been docked by AutoDockVina.

[0051] Executing step S100 includes the following specific steps:

[0052] NMDA targets were downloaded from the PubChem website, and 1 million ligand small molecule data sets to be screened were randomly selected from the REAL database.

[0053] Step S200: PLIP is used to calculate the binding pattern between the target to be screened and the small molecule, the binding pattern is encoded, and the distance matrix of the binding pattern is clustered using the DBSCAN algorithm to identify molecular groups with similar binding characteristics.

[0054] Combine Figure 2 , execute step S200, including the following specific steps:

[0055] Step S201: The target to be screened and the small molecule in step S100 are combined into a complex through PyMOL, hydrogen atoms are added to the small molecule, and excess water molecules in the complex are deleted.

[0056] Step S202: Calculate the binding modes of the complexes in step S201 using PLIP, which contains the following: Hydrogen Bonds, Hydrophobic Interactions, Salt Bridges, Pi-Cation Interactions, and Metal Complexes, a total of 5 non-covalent interaction types, and the corresponding residue types, residue numbers, and chain names.

[0057] Step S203: Encode the binding modes between the target and small molecules calculated in step S202, and use the Tanimoto coefficient to measure similarity, and calculate the distance matrix between the binding modes. The specific content includes:

[0058] First, align the Fasta sequences of the proteins to align the residue positions in the two protein chains. Then, compare whether there are the same non-covalent interaction types between the two protein-ligand complexes. After confirming the existence of the same non-covalent interaction types, further compare the residue pairs involved in these interactions. If the types, numbers, chain names, etc. of the two residue pairs are completely matched, then these interactions are considered the same. Next, digitally encode these binding modes, where the objects of encoding include:

[0059] Non-covalent interaction type: such as hydrogen bond, hydrophobic interaction, salt bridge, π-cation interaction, and metal complex.

[0060] Residue type: the type name of each residue.

[0061] Residue number: the number of the residue in the protein chain.

[0062] Chain name: the name of the protein chain to which the residue belongs.

[0063] The same type of non-covalent interaction under the same type of residue is encoded as 1.

[0064] Finally, measure the similarity between the two complexes by calculating the Tanimoto similarity, and further calculate the distance matrix. The formula for calculating the Tanimoto similarity is:

[0065]

[0066] Where A and B represent the encoded representations of the two complexes, ∩ represents the intersection, and ∪ represents the union.

[0067] The formula for calculating the distance matrix is:

[0068] Distance = 1 - Tanimoto Similarity

[0069] Step S204: Clustering the distance matrix of the binding mode in step S203 using the DBSCAN algorithm to identify groups of molecules with similar binding characteristics. The parameters of DBSCAN are:

[0070]

[0071] "eps":"0.5", which means the maximum distance of the neighborhood is 0.5.

[0072] "min_samples":"5", which means the minimum number of samples in a cluster is 5.

[0073] "metric":"precomputed", which means using the pre-computed distance matrix for calculation. DBSCAN will directly use the pre-computed distance matrix to determine the distance between points without calculating the distance itself.

[0074] Step S300: Selecting Top100 molecules for each class and calculating the hydrogen bond compensation number between proteins and ligands.

[0075] Performing step S300 includes the following specific steps:

[0076] Step S301: According to the AutoDockVina score, select Top100 molecules for each class in the clustering results in step S204, and use Gromacs command to generate a system file containing solvent. Process the target to be screened, small molecules and complexes respectively. Including the following specific content:

[0077] First, use PdbFixer to repair missing heavy atoms in the protein structure; second, use ACPYPE to generate small molecule topology files suitable for Gromacs; then, use Gromacs to generate topology files for the target protein to be screened and the complex, respectively, using the Amber03 force field. Finally, use Gromacs command to design a simulation box for the protein system and fill the box with water molecules as solvent, thus preparing for subsequent molecular dynamics simulation.

[0078] Step S302: Use PyMOL and the set threshold to calculate the number of hydrogen bonds in the solvated system file generated in step S301, and calculate the hydrogen bond compensation number.

[0079] Hydrogen bond compensation number = number of hydrogen bonds in LBD region complex - (number of hydrogen bonds in LBD region protein + number of hydrogen bonds in LBD region small molecule)

[0080] The tool used to find hydrogen bonds is the find_pairs function of Pymol, and the hyperparameters used are:

[0081]

[0082] "cutoff":"3.5", which represents the distance threshold of the hydrogen bond is 3.5.

[0083] "angle":"45", which represents the angle threshold of the hydrogen bond is 45°.

[0084] "mode":"1", which means that the angle condition will be considered when judging whether the hydrogen bond exists.

[0085] Step S400: Calculate the correlation between the residues involved in the binding mode of the Top100 molecules and the alpha carbon atoms of the key site.

[0086] Step S400 is executed, including the following specific steps:

[0087] Step S401: According to the AutoDockVina score, select the Top100 molecules of each class in the clustering results in step S204, use Gromacs to perform molecular dynamics simulation on the selected target, small molecule and complex, and generate topology file and trajectory file, including the following specific contents:

[0088] First, prepare the structure file and remove the water molecules. Second, generate the topology file, which uses Amber03 force field and specifies the water model as SPC / E. Next, add a box to the system and fill the solvent. Then, add ions to neutralize the system and ensure charge balance. After that, perform energy minimization to optimize the initial structure of the system. Next, perform NVT and NPT equilibration to adjust the temperature and pressure. Finally, perform the actual molecular dynamics simulation to complete the entire simulation process.

[0089] Step S402: Load the topology file and trajectory file generated in step S401 using MDAnalysis, calculate the Pearson correlation coefficient between the alpha carbon atoms of the key site and the relevant residues of the binding mode obtained in step S202, and then add these correlation coefficients to obtain the total correlation coefficient.

[0090] Step S500: Form an energy index according to the AutoDockVina score and the hydrogen bond compensation number, and form an activity index according to the correlation between the binding mode related residues and the alpha carbon atoms of the key site, and perform multi-objective frontier analysis on the selected TopK1 molecules to select the Pareto optimal solution set.

[0091] Step S500 is executed, including the following specific steps:

[0092] According to the hydrogen bond compensation number obtained in step S302 and the AutoDockVina score, the energy correlation index of each type of molecule is calculated by giving a weight of 0.2 and 0.8 respectively. According to the correlation value of the key site alpha carbon atom and the residue in the binding mode obtained in step S402, the activity correlation index is formed. In each type of molecule, based on the two dimensions of energy correlation and activity correlation, multi-objective frontier analysis is carried out, and the Pareto optimal solution set is selected as the final screening result.

[0093] It should be noted that the above content only illustrates the technical idea of the present application, and cannot limit the protection scope of the present application. For ordinary skilled persons in the art, without departing from the principles of the present application, a number of improvements and refinements can be made, which fall within the protection scope of the claims of the present application.

Claims

1. A virtual screening method for drug molecules based on binding mode analysis, characterized in that: The following steps are involved: S100: Prepare a target protein and a small molecule compound dataset to be screened, and use a molecular docking program to perform docking processing on the target protein and the small molecule compound dataset; S200: The docked target protein and small molecule compound datasets are merged into a complex using PyMOL. The binding mode of the complex is calculated using PLIP. The binding mode is encoded and the distance matrix is ​​calculated. The small molecule compound dataset is clustered using the distance matrix. S300: Selecting TopK1 molecules of each category in the clustering results of step S200, generating a system file containing a solvent for the target protein, small molecule compound, and complex to be screened after docking processing, and calculating the number of hydrogen bonds and hydrogen bond compensation number in the system file containing the solvent; S400: Select the TopK1 molecules of each category in the clustering results of step S200, perform molecular dynamics simulation on the target protein to be screened, small molecule compound and complex after docking processing, and calculate the correlation between the binding mode-related residues obtained in step S200 and the alpha carbon atom of the key site; S500: An energy correlation index is formed based on the hydrogen bond compensation number and AutoDockVina score, and an activity correlation index is formed based on the correlation between binding mode-related residues and alpha carbon atoms in key sites. Multi-objective analysis is performed on the selected TopK1 molecules to select the Pareto optimal solution set.

2. A virtual screening method for drug molecules based on binding mode analysis according to claim 1, characterized in that: In step S100, the molecular docking program is AutoDockVina.

3. The method for virtual screening of drug molecules based on binding mode analysis according to claim 1, characterized in that: In step S200 , the bonding modes include: Hydrogen Bonds, Hydrophobic Interactions, Salt Bridges, Pi-Cation Interactions, and Metal Complexes.

4. The method for virtual screening of drug molecules based on binding mode analysis according to claim 1, characterized in that: In the step S200 , the distance matrix of the binding patterns is calculated by using the Tanimoto coefficient to measure similarity and calculate the distance matrix between the binding patterns.

5. The method for virtual screening of drug molecules based on binding mode analysis according to claim 1, characterized in that: In step S200, the algorithm for clustering the distance matrix of the binding patterns is the DBSCAN algorithm.

6. The method for virtual screening of drug molecules based on binding mode analysis according to claim 1, characterized in that: In step S300, a system file containing a solvent is generated using a Gromacs command.

7. The method for virtual screening of drug molecules based on binding mode analysis according to claim 1, characterized in that: In step S300 , the number of hydrogen bonds and the number of hydrogen bond compensations in the system file containing the solvent are calculated using PyMOL and a set threshold.

8. The method for virtual screening of drug molecules based on binding mode analysis according to claim 1, characterized in that: described Step S400 specifically includes the following steps: Step S401: Based on the AutoDockVina score, select the TopK1 molecules of each category in the clustering results in step S200, perform molecular dynamics simulation on the docked target, small molecule and complex respectively using Gromacs, and generate topology files and trajectory files; Step S402: Use MDAnalysis to load the topology file and trajectory file generated in step S401, calculate the Pearson correlation coefficient between the alpha carbon atom of the key site and the binding mode-related residues obtained in step S200, and then add these correlation coefficients to obtain the total correlation coefficient.

9. A computer storage medium, characterized in that A computer program is stored, and when the computer program is executed by a processor, the drug molecule virtual screening method according to any one of claims 1 to 8 is implemented.

10. An electronic device, characterized in that: It comprises a memory and one or more processors, wherein the memory is used to store one or more programs; when the programs are executed by the one or more processors, the drug molecule virtual screening method according to any one of claims 1 to 8 is implemented.

Images