<jupyter_start><jupyter_text>Virtual ScreeningThis notebook aims to demonstrate how to use AutoDock Vina (via Smina) and Ledock to dock multiple molecules in the same protein target and binding site.  Content of this notebook1. Feching system and cleanup2. System Visualization 3. Docking with Smina    - Receptor preparation    - Ligand preparation    - Docking box definition    - Docking    - 3D visualization of docking results4. Docking with LeDock    - Receptor preparation    - Ligand preparation    - Docking box definition    - Docking    - DOK results file conversion to SDF    - 3D visualization of docking results<jupyter_code>from pymol import cmd
import py3Dmol

from openbabel import pybel

from rdkit import Chem
from rdkit.Chem import AllChem

import sys, os, random
sys.path.insert(1, 'utilities/')

from utils import getbox, generate_ledock_file, dok_to_sdf

import warnings
warnings.filterwarnings("ignore")

%config Completer.use_jedi = False
os.chdir('test/Virtual_Screening/')<jupyter_output><empty_output><jupyter_text>1. Feching system and cleanupImplementing Pymol is a simple way to download PDB structures. The user can launch this or any other Jupyter Dock's protocol by providing his or her own files.<jupyter_code>cmd.fetch(code='1X1R',type='pdb1')
cmd.select(name='Prot',selection='polymer.protein')
cmd.select(name='GDP',selection='organic')
cmd.save(filename='1X1R_clean.pdb',format='pdb',selection='Prot')
cmd.save(filename='1X1R_GDP.mol2',format='mol2',selection='GDP')
cmd.delete('all')<jupyter_output><empty_output><jupyter_text>2. System VisualizationA cool feature of Jupyter Dock is the posibility of visualize ligand-protein complexes and docking results into the notebbok. All this thanks to the powerful py3Dmol. Now the protein and ligand have been sanitized it would be recomended to vissualize the ligand-protein reference system.<jupyter_code>view = py3Dmol.view()
view.removeAllModels()
view.setViewStyle({'style':'outline','color':'black','width':0.1})

view.addModel(open('1X1R_clean.pdb','r').read(),format='pdb')
Prot=view.getModel()
Prot.setStyle({'cartoon':{'arrows':True, 'tubes':True, 'style':'oval', 'color':'white'}})
view.addSurface(py3Dmol.VDW,{'opacity':0.6,'color':'white'})

view.addModel(open('1X1R_GDP.mol2','r').read(),format='mol2')
ref_m = view.getModel()
ref_m.setStyle({},{'stick':{'colorscheme':'greenCarbon','radius':0.2}})

view.zoomTo()
view.show()<jupyter_output><empty_output><jupyter_text>3. Docking with SminaDespite the presence of Python bindings in AutoDock Vina 1.2.0, other tools that incorporate AutoDock Vina allow for cool features such as custom score functions (smina), fast execution (qvina), and the use of wider boxes (qvina-w). Jupyter Dock can run such binaries in a notebook, giving users more options.Smina is a fork of AutoDock Vina that is customized to better support scoring function development and high-performance energy minimization. Smina is maintained by <NAME> at the University of Pittsburgh and is not directly affiliated with the AutoDock project.>**Info:** The following cell contains an example of using Smina to run the current docking example. However, the executable files for qvina and qvina-w are available in the Jupyter Dock repo's bin directory. As a result, the user can use such a tool by adding the necessary cells or replacing the current docking engine. 3.1. Receptor preparationDespite the fact that Smina is a modified version of AutoDock Vina, the input file for a receptor in Smina can be either a PDBQT file or a PDB file with explicit hydrogens in all residues. At this point, we can make use ofca protein structure from Jupyter Dock's _**fix_protein()**_ function or implementing LePro (for more information, see the 2.1 section of the Molecular Docking notebook).<jupyter_code>!../../bin/lepro_linux_x86 {'1X1R_clean.pdb'}

os.rename('pro.pdb','1X1R_clean_H.pdb') # Output from lepro is pro.pdb, this line will change the name to '1X1R_clean_H.pdb'<jupyter_output><empty_output><jupyter_text>3.2. Ligand preparationThe ligand molecules in Virtual Screening protocols could come from a variety of sources (i.e. [ZINC15](https://zinc15.docking.org/), [PubChem](https://pubchem.ncbi.nlm.nih.gov/), [DrugBank](https://go.drugbank.com/), etc) and diverse formats (i.e. SDF, PDB, MOL, MOL2, SMILES, etc). The cell below depicts one of the simplest approaches to using molecules from SMILES codes. Users, however, can use any known chemical format in their molecules thanks to the use of PyBel and RDKit.Regardless of the format of the molecules or the differences between docking algorithms, ligand preparation must achieve at least the following objectives:- Set the proper protonation and tautomeric state to the molecules.- Provide a valid 3D structure to initialize the conformational search.Jupyter Dock generates protonated neutral molecules and performs energy minimization under the MMFF94s force field using the pybel functions _make3D_ and _localopt_, respectively.It is recommended that users tune the pybel settings or use other preparation methods before running the molecular docking for special ligand needs.<jupyter_code>smiles=['C1=NC(=C2C(=N1)N(C=N2)CCOCP(=O)(O)O)N',
        'C1=NC(=C2C(=N1)N(C=N2)[C@H]3[C@@H]([C@@H]([C@H](O3)COP(=O)(O)O)O)O)N',
        'C[C@@H](C1=CC2=C(C(=CC=C2)Cl)C(=O)N1C3=CC=CC=C3)NC4=NC=NC5=C4NC=N5',
        'C1=NC(=C2C(=N1)N(C=N2)C3C(C(C(O3)CO)O)O)N',
        'C[C@H](CN1C=NC2=C(N=CN=C21)N)OC[P@@](=O)(N[C@@H]',
        'C1=NC(=C2C(=N1)N(C=N2)[C@H]3[C@@H]([C@@H]([C@H](O3)C[C@H](CC[C@@H](C(=O)O)N)N)O)O)N',
        'C[S+](CC[C@@H](C(=O)O)N)C[C@@H]1[C@H]([C@H]([C@@H](O1)N2C=NC3=C(N=CN=C32)N)O)O',
        'CN1C=C(C(=N1)OC)NC2=C3C(=NC(=N2)N4C[C@H]([C@@H](C4)F)NC(=O)C=C)N(C=N3)C',
        'C1COC[C@@H]1NC2=C3C(=NC=N2)N(C=N3)[C@H]4[C@@H]([C@@H]([C@H](O4)CO)O)O']
out=pybel.Outputfile(filename='InputMols.mol2',format='mol2',overwrite=True)
for index,smi in enumerate(smiles):
    mol=pybel.readstring(string=smi,format='smiles')
    mol.title='mol_'+str(index)
    mol.make3D('mmff94s')
    mol.localopt(forcefield='mmff94s', steps=500)
    out.write(mol)
out.close()<jupyter_output><empty_output><jupyter_text>3.3. Docking box definitionIn this example, the natural ligand GDP of 1X1R will be used as a reference for box definition (Section 4.3 of the Molecular Docking Notebook contains more information about the _**getbox()**_ function).<jupyter_code>cmd.load(filename='1X1R_clean_H.pdb',format='pdb',object='prot') #Not needed but as reference of the system
cmd.load(filename='1X1R_GDP.mol2',format='mol2',object='lig')

center,size=getbox(selection='lig',extending=6.0,software='vina')
cmd.delete('all')
print(center)
print(size)<jupyter_output>{'center_x': 3.4204999953508377, 'center_y': 9.91599988937378, 'center_z': 11.27299976348877}
{'size_x': 19.56700000166893, 'size_y': 18.30399990081787, 'size_z': 23.20599937438965}<jupyter_text>3.4. DockingJupyter Dock comes with the smina executable for Linux and Mac OS. By running the binary file, the parameters can be accessed.>**Hint:** Smina also has the ability to use a multimodel file as input and then save the results to a single multimodel output file. Furthermore, the input and output molecule files can be tar.gz. sdf files, which is very useful when working with large databases.<jupyter_code>!../../bin/smina -r {'1X1R_clean_H.pdb'} -l {'InputMols.mol2'} -o {'1X1R_lig_smina_out.sdf'} --center_x {center['center_x']} --center_y {center['center_y']} --center_z {center['center_z']} --size_x {size['size_x']} --size_y {size['size_y']} --size_z {size['size_z']} --exhaustiveness 8 --num_modes 5<jupyter_output>_______  _______ _________ _        _______ 
  (  ____ \(       )\__   __/( (    /|(  ___  )
  | (    \/| () () |   ) (   |  \  ( || (   ) |
  | (_____ | || || |   | |   |   \ | || (___) |
  (_____  )| |(_)| |   | |   | (\ \) ||  ___  |
        ) || |   | |   | |   | | \   || (   ) |
  /\____) || )   ( |___) (___| )  \  || )   ( |
  \_______)|/     \|\_______/|/    )_)|/     \|


smina is based off AutoDock Vina. Please cite appropriately.

Weights      Terms
-0.035579    gauss(o=0,_w=0.5,_c=8)
-0.005156    gauss(o=3,_w=2,_c=8)
0.840245     repulsion(o=0,_c=8)
-0.035069    hydrophobic(g=0.5,_b=1.5,_c=8)
-0.587439    non_dir_h_bond(g=-0.7,_b=0,_c=8)
1.923        num_tors_div

Using random seed: 1174223054

0%   10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
***************************************************

mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
[...]<jupyter_text>3.5. 3D visualization of docking resultsAs with the system visualization (section 2), the docking results can be inspected and compared to the reference structure (if one exists). Smina saves the "minimizedAffinity" information corresponding to the docking score as the molecule's attribute.<jupyter_code>view = py3Dmol.view()
view.removeAllModels()
view.setViewStyle({'style':'outline','color':'black','width':0.1})

view.addModel(open('1X1R_clean_H.pdb','r').read(),'pdb')
Prot=view.getModel()
Prot.setStyle({'cartoon':{'arrows':True, 'tubes':True, 'style':'oval', 'color':'white'}})
view.addSurface(py3Dmol.VDW,{'opacity':0.8,'color':'white'})


view.addModel(open('1X1R_GDP.mol2','r').read(),'mol2')
ref_m = view.getModel()
ref_m.setStyle({},{'stick':{'colorscheme':'greenCarbon','radius':0.2}})


poses=Chem.SDMolSupplier('1X1R_lig_smina_out.sdf',True)
for p in list(poses)[::5]:
    pose_1=Chem.MolToMolBlock(p)
    print(p.GetProp('_Name'),'Score: {}'.format(p.GetProp('minimizedAffinity')))
    color = ["#"+''.join([random.choice('0123456789ABCDEF') for j in range(6)])]
    view.addModel(pose_1,'mol')
    z= view.getModel()
    z.setStyle({},{'stick':{'color':color[0],'radius':0.05,'opacity':0.6}})

view.zoomTo()
view.show()<jupyter_output>mol_0 Score: -7.07914
mol_1 Score: -9.79016
mol_2 Score: -9.37383
mol_3 Score: -8.05686
mol_4 Score: -7.05481
mol_5 Score: -10.09781
mol_6 Score: -9.98733
mol_7 Score: -8.02271
mol_8 Score: -8.91059<jupyter_text>4. Docking with LeDockLeDock is designed for fast and accurate flexible docking of small molecules into a protein. It achieves a pose-prediction accuracy of greater than 90% on the Astex diversity set and takes about 3 seconds per run for a drug-like molecule. It has led to the discovery of novel kinase inhibitors and bromodomain antagonists from high-throughput virtual screening campaigns. It directly uses SYBYL Mol2 format as input for small molecules. 4.1. Receptor preparationIn LeDock, the input file for a receptor is a PDB file with explicit hydrogens in all residues. LePro was created as a tool for preparing protein structures for docking with LeDock. Thus, at this stage, we can use the file ater sanitization steps as well as a protein structure from Jupyter Dock's _**fix_protein()**_ function or implementing LePro (for more information, see the 2.1 section of the Molecular Docking notebook). 4.2. Ligand preparationMOL2 is the ligand input format in LeDock. LeDock, on the other hand, cannot accept a multimodel MOL2 file as an input. To use LeDock, we can run the preparation in the same way as Smina (section 3.2), but generate one MOL2 file per ligand.<jupyter_code>smiles=['C1=NC(=C2C(=N1)N(C=N2)CCOCP(=O)(O)O)N',
        'C1=NC(=C2C(=N1)N(C=N2)[C@H]3[C@@H]([C@@H]([C@H](O3)COP(=O)(O)O)O)O)N',
        'C[C@@H](C1=CC2=C(C(=CC=C2)Cl)C(=O)N1C3=CC=CC=C3)NC4=NC=NC5=C4NC=N5',
        'C1=NC(=C2C(=N1)N(C=N2)C3C(C(C(O3)CO)O)O)N',
        'C[C@H](CN1C=NC2=C(N=CN=C21)N)OC[P@@](=O)(N[C@@H]',
        'C1=NC(=C2C(=N1)N(C=N2)[C@H]3[C@@H]([C@@H]([C@H](O3)C[C@H](CC[C@@H](C(=O)O)N)N)O)O)N',
        'C[S+](CC[C@@H](C(=O)O)N)C[C@@H]1[C@H]([C@H]([C@@H](O1)N2C=NC3=C(N=CN=C32)N)O)O',
        'CN1C=C(C(=N1)OC)NC2=C3C(=NC(=N2)N4C[C@H]([C@@H](C4)F)NC(=O)C=C)N(C=N3)C',
        'C1COC[C@@H]1NC2=C3C(=NC=N2)N(C=N3)[C@H]4[C@@H]([C@@H]([C@H](O4)CO)O)O']
for index,smi in enumerate(smiles):
    mol=pybel.readstring(string=smi,format='smiles')
    mol.title='mol_'+str(index)
    mol.make3D('mmff94s')
    mol.localopt(forcefield='mmff94s', steps=500)
    out=pybel.Outputfile(filename='ledock_inputfiles/'+'mol_'+str(index)+'.mol2',format='mol2',overwrite=True)
    out.write(mol)
    out.close()<jupyter_output><empty_output><jupyter_text>4.3. Docking box definitionThis step can be completed in the same manner as the Smina box definition (section 3.3). To obtain the identical box from Smina docking but in LeDock format, the user only needs to change the parameter "software" from "vina" to "ledock." >**Info:** The implementation of the _**getbox()**_ function allows for the easy replication of binding sites between AutoDock Vina and LeDock, with the goal of replicating and comparing results between both programs.<jupyter_code>cmd.load(filename='1X1R_clean_H.pdb',format='pdb',object='prot')
cmd.load(filename='1X1R_GDP.mol2',format='mol2',object='lig')

X,Y,Z=getbox(selection='lig',extending=6.0,software='ledock')
cmd.delete('all')
print(X)
print(Y)
print(Z)<jupyter_output>{'minX': -6.363000005483627, 'maxX': 13.203999996185303}
{'minY': 0.7639999389648438, 'maxY': 19.067999839782715}
{'minZ': -0.3299999237060547, 'maxZ': 22.875999450683594}<jupyter_text>4.4 DockingAside from one file per ligand, Ledock requires the user to create a file containing the docking paths for the ligands (ligand.list in Jupyter Dock). This file can be easily created by using the following cell and the _**generate_ledock_file()**_ function.<jupyter_code>l_list=[]
for file in os.listdir('ledock_inputfiles/'):
    if 'mol2' in file:
        l_list.append('ledock_inputfiles/'+file+'\n')
l_list
generate_ledock_file(receptor='1X1R_clean_H.pdb',l_list=l_list,
                     l_list_outfile='ligand.list',
                     x=[X['minX'],X['maxX']],
                     y=[Y['minY'],Y['maxY']],
                     z=[Z['minZ'],Z['maxZ']],
                     n_poses=10,
                     rmsd=1.0,
                     out='dock.in')<jupyter_output><empty_output><jupyter_text>Once all of the docking parameters have been entered into a configuration file (dock.in), everything is ready to go.<jupyter_code>!../../bin/ledock_linux_x86 {'dock.in'}<jupyter_output><empty_output><jupyter_text>4.5. DOK results file conversion to SDFLeDock produces a file with the dok extension that contains docking properties in the same way that a pdb file does. Despite this, the dok file is not widely used for representing chemical structures. As a result, Jupyter Dock is capable of converting dok files to the widely used sdf format. Jupyter Dock will save the "Pose" and "Score" results as molecule attributes while preserving the chemical features (more information in section 6.5 of the Molecular Docking notebook).<jupyter_code>for file in os.listdir('ledock_inputfiles/'):
    if '.dok' in file:
        os.rename('ledock_inputfiles/'+file,'ledock_outfiles/'+file)
        dok_to_sdf(dok_file='ledock_outfiles/'+file,output='ledock_outfiles/'+file.replace('dok','sdf'))<jupyter_output><empty_output><jupyter_text>4.5. 3D visualization of docking resultsAs with the system visualization (section 2), the docking results can be inspected and compared to the reference structure (if one exists).<jupyter_code>view = py3Dmol.view()
view.removeAllModels()
view.setViewStyle({'style':'outline','color':'black','width':0.1})

view.addModel(open('1X1R_clean_H.pdb','r').read(),'pdb')
Prot=view.getModel()
Prot.setStyle({'cartoon':{'arrows':True, 'tubes':True, 'style':'oval', 'color':'white'}})
view.addSurface(py3Dmol.VDW,{'opacity':0.8,'color':'white'})


view.addModel(open('1X1R_GDP.mol2','r').read(),'mol2')
ref_m = view.getModel()
ref_m.setStyle({},{'stick':{'colorscheme':'greenCarbon','radius':0.2}})

for file in os.listdir('ledock_outfiles/'):
    if 'sdf' in file:
        pose_1=Chem.SDMolSupplier('ledock_outfiles/'+file,False)[0]
        p=Chem.MolToMolBlock(pose_1)
        print('Name: {} | Pose: {} | Score: {}'.format(file.split('.')[0],pose_1.GetProp('Pose'),pose_1.GetProp('Score')))
        color = ["#"+''.join([random.choice('0123456789ABCDEF') for j in range(6)])]
        view.addModel(p,'mol')
        z= view.getModel()
        z.setStyle({},{'stick':{'color':color[0],'radius':0.05,'opacity':0.6}})

view.zoomTo()
view.show()<jupyter_output>Name: mol_0 | Pose: 2 | Score: -7.91
Name: mol_1 | Pose: 5 | Score: -9.55
Name: mol_2 | Pose: 3 | Score: -7.84
Name: mol_3 | Pose: 1 | Score: -7.08
Name: mol_4 | Pose: 1 | Score: -7.51
Name: mol_5 | Pose: 2 | Score: -9.46
Name: mol_6 | Pose: 4 | Score: -9.99
Name: mol_7 | Pose: 1 | Score: -8.42
Name: mol_8 | Pose: 2 | Score: -7.16
{'fraction_non_alphanumeric': 0.16085838962281068, 'fraction_numerical': 0.04417478629924084, 'mean_word_length': 5.007181328545781, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 56, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '124733', 'license': 'MIT', 'n_tokens_mistral': 6654, 'n_tokens_neox': 6100, 'n_words': 1489, 'path': '.ipynb_checkpoints/2.-Virtual_Screening-checkpoint.ipynb', 'repo_name': 'amoliu/Jupyter_Dock', 'size': '16733'}