Drug molecule virtual screening method based on binding mode analysis

Abstract

The invention discloses a drug molecule virtual screening method based on combination mode analysis, and belongs to the field of artificial intelligence drug discovery. Using PLIP to calculate binding modes between the to-be-screened target and the small molecules, performing coding, using a DBSCAN algorithm to cluster the binding modes, and identifying compound populations of similar binding modes; in the screening stage, two key indexes are introduced: an energy index is formed according to a hydrogen bond compensation number and an AutoDock Vina score, the correlation between related residues and key site alpha carbon atoms in a binding mode is used as an activity index for evaluating the binding stability and characteristic activity between a target and a small molecule, and multi-target frontier analysis is performed on candidate molecules; according to the method, the efficiency and precision of virtual screening are remarkably improved, detailed analysis on intermolecular interaction is enhanced, the method is suitable for small-scale and large-scale molecular screening, the drug research and development cost is reduced, the drug research and development time is shortened, the diversity of virtual screening results is enriched, and potential candidate molecules are provided for subsequent experiments.

Description

Technical Field

[0001] The present invention belongs to the technical field of artificial intelligence drug discovery, and specifically relates to a virtual screening method for drug molecules based on binding pattern analysis. Background Art

[0002] Since its birth, virtual screening technology has become an indispensable part of modern drug discovery. Through computer simulation, scientists can screen a large number of small molecule compounds in a short period of time and predict their binding affinity with the target protein, thereby providing 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 significantly reduce screening time, reduce costs, and improve screening efficiency.

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

[0004] The introduction of binding mode analysis has brought a qualitative leap in 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 analytical 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, thereby improving the activity and selectivity of drugs. 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 set of candidates for subsequent experimental screening.

[0005] Although binding mode analysis provides more depth and precision for virtual screening, it is still not enough to fully optimize the diversity and accuracy of screening results by relying on this technology alone. For example, existing methods may not be able to fully evaluate the combined effects 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 shortcomings in order to achieve more comprehensive and accurate screening results. Summary of the invention

[0006] To solve the above problems, the present invention discloses a virtual screening method for drug molecules based on binding mode analysis. PLIP is used to calculate the binding mode between the target to be screened and the small molecule, the binding mode is encoded, and the binding mode is clustered using the DBSCAN algorithm to identify groups of compounds with similar binding modes. In the screening stage, the present invention further introduces two key indicators: an energy index is formed based on the hydrogen bond compensation number and the AutoDockVina score, and the correlation between the relevant residues in the binding mode and the alpha carbon atom at the key site is used as an activity indicator to evaluate the binding stability between the target and the small molecule and the characteristic activity of the binding mode, and to perform multi-target frontier analysis on the candidate molecules.

[0007] To achieve the above object, the technical solution of the present invention is as follows:

[0008] The present invention provides a virtual screening method for drug molecules based on binding mode analysis, comprising the following steps:

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

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

[0011] Step S300: Select TopK1 molecules in each category and calculate the hydrogen bond compensation number between the protein and the ligand;

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

[0013] Step S500: An energy index is formed according to the hydrogen bond compensation number and the AutoDockVina score, and an activity index is formed according to the correlation between the binding mode-related residues and the alpha carbon atom of the key site, and a multi-objective frontier analysis is performed on the selected TopK1 molecules to select the Pareto optimal solution set.

[0014] Furthermore, step S100 specifically includes the following steps:

[0015] Select a target protein to be screened and prepare the relevant small molecule compound data set; use AutoDockVina to dock the target with the small molecule data set and generate a conformation file after docking.

[0016] Furthermore, step S200 specifically includes the following steps:

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

[0018] Step S202: Use PLIP to calculate the binding mode of the complex in step S201. The binding mode includes five non-covalent interactions, namely Hydrogen Bonds, Hydrophobic Interactions, Salt Bridges, Pi-Cation Interactions and Metal Complexes, as well as the corresponding residue types, residue numbers and chain names.

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

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

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

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

[0023] Hydrophobic Interactions: Hydrophobic interactions are non-covalent interactions caused by hydrophobic molecules avoiding water and clustering in space.

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

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

[0026] 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.

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

[0028] Step S301: According to the AutoDockVina score, select the TopK1 molecules of each category in the clustering results in step S204, use the Gromacs command to generate a system file containing the solvent, and process the target to be screened, small molecules and complexes respectively.

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

[0030] Furthermore, step S400 specifically includes the following steps:

[0031] Step S401: According to the AutoDockVina score, select the TopK1 molecules of each category in the clustering results in step S204, use Gromacs to perform molecular dynamics simulation on the docked target, small molecule and complex, and generate a topology file and a trajectory file;

[0032] 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 S202, and then add these correlation coefficients to obtain the total correlation coefficient.

[0033] Furthermore, step S500 specifically includes the following steps:

[0034] Energy indices were formed based on the hydrogen bond compensation number and the AutoDockVina score, and activity indices were formed based on the correlation between binding mode-related residues and alpha carbon atoms at key sites. Multi-objective frontier analysis was performed on the selected TopK1 molecules to select the Pareto optimal solution set.

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

[0036] 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. First, a target to be screened and a small molecule data set that 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 molecule, and the binding mode is encoded. On this basis, the distance matrix of the binding mode is clustered using the DBSCAN algorithm to identify molecular groups with similar binding characteristics; then, in the preliminary molecular screening stage, TopK1 molecules are selected from each cluster according to the AutoDockVina score, and the hydrogen bond compensation number between the target to be screened and the small molecule is further calculated, and the correlation between the residues involved in the binding mode of the TopK1 molecule and the alpha carbon atom of the key site is evaluated; finally, an energy index is formed according to the hydrogen bond compensation number and the AutoDockVina score, and the correlation between the relevant residues in the binding mode and the alpha carbon atom of the key site is used as the activity index, and the selected TopK1 molecules are subjected to multi-target frontier analysis, and the Pareto optimal solution set is selected as the final screening result.

[0037] The method of the present invention 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 correlation between the hydrogen bond compensation number and the alpha carbon atom at the key site, the analysis of intermolecular interactions is further optimized. Ultimately, this systematic and refined analysis method, combined with multi-target frontier analysis, not only significantly improves the effect of virtual screening, and is suitable for small-scale and large-scale molecular screening, but also reduces the cost and time of drug development, and enriches the diversity of virtual screening results, achieving more comprehensive and accurate screening results, and providing more potential candidate molecules for subsequent experiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 The figure is an overall flow chart of the virtual screening method of drug molecules based on binding mode analysis of the present invention;

[0039] Figure 2 Computational flow chart for binding mode analysis. DETAILED DESCRIPTION

[0040] The present invention will be further explained below in conjunction with the accompanying drawings and specific embodiments. It should be understood that the following specific embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention.

[0041] The NMDA target was chosen as an exemplary target for this method.

[0042] Non-covalent interactions are explained as follows:

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

[0044] Hydrophobic Interactions: Hydrophobic interactions are non-covalent interactions caused by hydrophobic molecules avoiding water and clustering 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 attractions 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] Combination 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 the 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 mode between the target to be screened and the small molecule, the binding mode is encoded, and the distance matrix of the binding mode is clustered using the DBSCAN algorithm to identify molecular groups with similar binding characteristics.

[0054] Combination 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 redundant water molecules in the complex are deleted.

[0056] Step S202: Use PLIP to calculate the binding mode of the complex in step S201. The binding mode includes five non-covalent interactions, namely Hydrogen Bonds, Hydrophobic Interactions, Salt Bridges, Pi-Cation Interactions and Metal Complexes, as well as the corresponding residue types, residue numbers and chain names.

[0057] Step S203: Encode the binding mode between the target and the small molecule calculated in step S202, and use the Tanimoto coefficient to measure the similarity and calculate the distance matrix between the binding modes. The specific contents include the following:

[0058] First, the Fasta sequences of the proteins are aligned to align the residue positions in the two protein chains. Next, the two protein-ligand complexes are compared to see if they have the same type of non-covalent interactions. After confirming the presence of the same type of non-covalent interactions, the residue pairs involved in these interactions are further compared. If the type, number, chain name, etc. of the two residue pairs match exactly, then these interactions are considered to be the same. Next, these binding patterns are digitally encoded, where the objects encoded include:

[0059] Types of non-covalent interactions: such as hydrogen bonds, hydrophobic interactions, salt bridges, π-cation interactions, and metal complexes.

[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 residues in the same type of non-covalent interaction are coded as 1.

[0064] Finally, the similarity between the two complexes is measured by calculating the Tanimoto similarity, and the distance matrix is ​​further calculated. The calculation formula of Tanimoto similarity is:

[0065]

[0066] Among them, A and B represent the encoding representations of two complexes respectively, ∩ represents the intersection, and ∪ represents the union.

[0067] The calculation formula of the distance matrix is:

[0068] Distance=1-Tanimoto Similarity

[0069] Step S204: clustering the distance matrix of the binding patterns 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", 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", means using a precomputed distance matrix for calculation. DBSCAN will directly use the precomputed distance matrix to determine the distance between points without calculating the distance itself.

[0074] Step S300: Select the top 100 molecules in each category and calculate the hydrogen bond compensation number between the protein and the ligand.

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

[0076] Step S301: According to the AutoDockVina score, select the top 100 molecules of each category in the clustering results in step S204, use the Gromacs command to generate a system file containing the solvent, and process the target to be screened, small molecules and complexes respectively. The specific contents include the following:

[0077] First, PdbFixer was used to repair the missing heavy atoms in the protein structure. Second, ACPYPE was used to generate a small molecule topology file suitable for Gromacs. Then, Gromacs was used to generate the topology files of the target protein and complex to be screened, respectively, using the Amber03 force field. Finally, the Gromacs command was used to design a simulation box for the protein system, and water molecules were used to fill the box as a solvent to prepare for the subsequent molecular dynamics simulation.

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

[0079] Hydrogen bond compensation number = number of hydrogen bonds in the LBD region complex - (number of hydrogen bonds in the LBD region protein + number of hydrogen bonds in the 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 means the distance threshold of hydrogen bond is 3.5.

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

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

[0085] Step S400: Calculate the correlation between the residues involved in the Top 100 molecular binding modes and the alpha carbon atoms at the key sites.

[0086] Executing step S400 includes the following specific steps:

[0087] Step S401: According to the AutoDockVina score, select the top 100 molecules of each category in the clustering results in step S204, use Gromacs to perform molecular dynamics simulation on the selected targets to be screened, small molecules and complexes, and generate topology files and trajectory files, including the following specific contents:

[0088] First, prepare the structure file and remove the water molecules. Next, generate the topology file, using the Amber03 force field and specifying the water model as SPC / E. Next, add boxes to the system and fill it with 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 balancing to adjust the temperature and pressure. Finally, perform the actual molecular dynamics simulation to complete the entire simulation process.

[0089] 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 S202, and then add these correlation coefficients to obtain the total correlation coefficient.

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

[0091] Executing step S500 includes the following specific steps:

[0092] According to the hydrogen bond compensation number and AutoDockVina score obtained in step S302, the weights of 0.2 and 0.8 are assigned respectively to calculate the energy correlation index of each type of molecule. According to the correlation value between 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, a multi-objective frontier analysis is performed to select the Pareto optimal solution set as the final screening result.

[0093] It should be noted that the above content only illustrates the technical idea of ​​the present invention and cannot be used to limit the protection scope of the present invention. For ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications all fall within the protection scope of the claims of the present invention.

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 data set to be screened, and use a molecular docking program to perform docking processing on the target protein and the small molecule compound data set; 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, and the small molecule compound dataset is clustered using the distance matrix; S300: Selecting TopK1 molecules of each category in the clustering results in 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 the number of hydrogen bond compensations in the system file containing the solvent; S400: Select the TopK1 molecules of each category in the clustering results in 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 the AutoDockVina score, and an activity correlation index is formed based on the correlation between the binding mode-related residues and the alpha carbon atoms in the key sites. A 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 the step S100, the molecular docking program is AutoDockVina.

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

4. A virtual screening method for 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 the 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. A virtual screening method for drug molecules based on binding mode analysis according to claim 1, characterized in that: In step S300, a system file containing solvent is generated using a Gromacs command.

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

8. The method for virtual screening of drug molecules based on binding mode analysis according to claim 1, characterized in that: Said Step S400 specifically includes the following steps: Step S401: According to the AutoDockVina score, select the TopK1 molecules of each category in the clustering results in step S200, use Gromacs to perform molecular dynamics simulation on the docked target, small molecule and complex, and generate a topology file and a trajectory file; 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 method for virtual screening of drug molecules as claimed in 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 as described in any one of claims 1 to 8 is implemented.

Images