Meeko reads an RDKit molecule object and writes a PDBQT string (or file) for AutoDock-Vina and AutoDock-GPU. It converts the docking output to RDKit molecules and SD files, without loss of bond orders.
Meeko is developed by the Forli lab at the Center for Computational Structural Biology (CCSB) at Scripps Research.
Meeko does not calculate 3D coordinates or assign protonation states. Input molecules must have explicit hydrogens.
Sampling of macrocycle conformers (paper) is enabled by default.
SDF format strongly preferred over Mol2. See this discussion, and this one, also this, and this. and RDKit issues 1755 and 917. So, what could go wrong? For example, reading Mol2 files from ZINC led to incorrect net charge of some molecules.
Class MoleculePreparation no longer has method write_pdbqt_string().
Instead, MoleculePreparation.prepare() returns a list of MoleculeSetup objects
that must be passed, individually, to PDBQTWriterLegacy.write_string().
from meeko import MoleculePreparation
from meeko import PDBQTWriterLegacy
preparator = MoleculePreparation()
mol_setups = preparator.prepare(rdkit_molecule_3D_with_Hs)
for setup in mol_setups:
pdbqt_string, is_ok, error_msg = PDBQTWriterLegacy.write_string(setup)
if is_ok:
print(pdbqt_string, end="")Argument keep_nonpolar_hydrogens is replaced by merge_these_atom_types, both in the Python
interface and for script mk_prepare_ligand.py.
The default is merge_these_atom_types=("H",), which
merges hydrogens typed "H", keeping the current default behavior.
To keep all hydrogens, set merge_these_atom_types to an empty
list when initializing MoleculePreparation, or pass no atom types
to --merge_these_atom_types from the command line:
mk_prepare_ligand.py -i molecule.sdf --merge_these_atom_types- Python (>=3.5)
- Numpy
- Scipy
- RDKit
- ProDy (optionally, for covalent docking)
Conda or Miniconda can install the dependencies:
conda install -c conda-forge numpy scipy rdkit
pip install prody # optional. pip recommended at http://prody.csb.pitt.edu/downloads/$ pip install meekoIf using conda, pip installs the package in the active environment.
You'll get the develop branch, which may be ahead of the latest release.
$ git clone https://github.com/forlilab/Meeko
$ cd Meeko
$ pip install .Optionally include --editable. Changes in the original package location
take effect immediately without the need to run pip install . again.
$ pip install --editable .AutoDock-GPU and Vina read molecules in the PDBQT format. These can be prepared by Meeko from SD files, or from Mol2 files, but SDF is strongly preferred.
mk_prepare_ligand.py -i molecule.sdf -o molecule.pdbqt
mk_prepare_ligand.py -i multi_mol.sdf --multimol_outdir folder_for_pdbqt_filesAutoDock-GPU and Vina write docking results in the PDBQT format. The DLG output from AutoDock-GPU contains docked poses in PDBQT blocks. Meeko generates RDKit molecules from PDBQT files (or strings) using the SMILES string in the REMARK lines. The REMARK lines also have the mapping of atom indices between SMILES and PDBQT. SD files with docked coordinates are written from RDKit molecules.
mk_export.py molecule.pdbqt -o molecule.sdf
mk_export.py vina_results.pdbqt -o vina_results.sdf
mk_export.py autodock-gpu_results.dlg -o autodock-gpu_results.sdfMaking RDKit molecules from SMILES is safer than guessing bond orders from the coordinates, specially because the PDBQT lacks hydrogens bonded to carbon. As an example, consider the following conversion, in which OpenBabel adds an extra double bond, not because it has a bad algorithm, but because this is a nearly impossible task.
$ obabel -:"C1C=CCO1" -o pdbqt --gen3d | obabel -i pdbqt -o smi
[C]1=[C][C]=[C]O1from meeko import MoleculePreparation
from meeko import PDBQTWriterLegacy
from rdkit import Chem
input_molecule_file = "example/BACE_macrocycle/BACE_4.sdf"
# there is one molecule in this SD file, this loop iterates just once
for mol in Chem.SDMolSupplier(input_molecule_file, removeHs=False):
preparator = MoleculePreparation()
mol_setups = preparator.prepare(mol)
for setup in mol_setups:
setup.show() # optional
pdbqt_string = PDBQTWriterLegacy.write_string(setup)At this point, pdbqt_string can be written to a file for
docking with AutoDock-GPU or Vina, or passed directly to Vina within Python
using set_ligand_from_string(pdbqt_string). For context, see
the docs on using Vina from Python.
from meeko import PDBQTMolecule
from meeko import RDKitMolCreate
fn = "autodock-gpu_results.dlg"
pdbqt_mol = PDBQTMolecule.from_file(fn, is_dlg=True, skip_typing=True)
rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)The length of rdkitmol_list is one if there are no sidechains and only one
ligand was docked.
If multiple ligands and/or sidechains are docked simultaneously, each will be
an individual RDKit molecule in rdkitmol_list.
Sidechains are truncated at the C-alpha.
Note that docking multiple
ligands simultaneously is only available in Vina, and it differs from docking
multiple ligands one after the other. Each failed creation of an RDKit molecule
for a ligand or sidechain results in a None in rdkitmol_list.
For Vina's output PDBQT files, omit is_dlg=True.
pdbqt_mol = PDBQTMolecule.from_file("vina_results.pdbqt", skip_typing=True)
rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)When using Vina from Python, the output string can be passed directly.
See the docs
for context on the v object.
vina_output_string = v.poses()
pdbqt_mol = PDBQTMolecule(vina_output_string, is_dlg=True, skip_typing=True)
rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)This is useful for saving and loading configuration files with json.
import json
from meeko import MoleculePreparation
mk_config = {"hydrate": True} # any arguments of MoleculePreparation.__init__
print(json.dumps(mk_config), file=open('mk_config.json', 'w'))
with open('mk_config.json') as f:
mk_config = json.load(f)
preparator = MoleculePreparation.from_config(mk_config)The command line tool mk_prepare_ligand.py can read the json files using
option -c or --config.
Here we create an instance of MoleculePreparation that attaches pseudo waters to the ligand (see paper on hydrated docking), keeps macrocycles rigid (generally a bad idea), and keeps conjugated bonds and amide bonds rigid. By default, most amides are kept rigid, except for tertiary amides with different substituents on the nitrogen.
preparator = MoleculePreparation(
hydrate=True,
rigid_macrocycles=True,
rigidify_bonds_smarts = ["C=CC=C", "[CX3](=O)[NX3]"],
rigidify_bonds_indices = [(1, 2), (0, 2)],
)The same can be done with the command line script. Note that indices of the atoms in the SMARTS are 0-based for the Python API but 1-based for the command line script:
mk_prepare_ligand.py\
--hydrate\
--rigid_macrocycles\
--rigidify_bonds_smarts "C=CC=C"\
--rigidify_bonds_indices 2 3\
--rigidify_bonds_smarts "[CX3](=O)[NX3]"\
--rigidify_bonds_indices 1 3The input ligand must be the product of the reaction and contain the
atoms of the flexible sidechain up to (and including) the C-alpha.
For example, if the ligand is an acrylamide (smiles: C=CC(=O)N) reacting
with a cysteine (sidechain smiles: CCS), then the input ligand for
Meeko must correspond to smiles CCSCCC(=O)N.
Meeko will align the ligand atoms that match the C-alpha and C-beta of
the protein sidechain. Options --tether_smarts and --tether_smarts_indices
define these atoms. For a cysteine, --tether_smarts "SCC" and
--tether_smarts_indices 3 2 would work, although it is safer to define
a more spefic SMARTS to avoid matching the ligand more than once. The first
index (3 in this example) defines the C-alpha, and the second index defines
the C-beta.
For the example in this repository, which is based on PDB entry 7aeh, the following options prepare the covalently bound ligand for tethered docking:
cd example/covalent_docking
mk_prepare_ligand.py\
-i ligand_including_cys_sidechain.sdf\
--receptor protein.pdb\
--rec_residue ":CYS:6"\
--tether_smarts "NC(=O)C(O)(C)SCC"\
--tether_smarts_indices 9 8\
-o prepared.pdbqtFollow wk_prepare_receptor.py instructions and run with --pdb.
The goal of this step is to perform essential fixes to the protein
(such as missing atoms), to add hydrogens, and to follow the Amber
naming scheme for atoms and residues, e.g., HIE or HID
instead of HIS.
Here, wk.pdb was written by waterkit.
Here, wk.pdb was written by waterkit. The example below will center a gridbox of specified size on the given reactive residue.
$ mk_prepare_receptor.py\
--pdb wk.pdb\
-o receptor.pdbqt\
--flexres " :ARG:348"\
--reactive_flexres " :SER:308"
--reactive_flexres " :SER:308"\
--box_center_on_reactive_res\
--box_size 40 40 40 # x y z (angstroms)A manual box center can be specified with --box_center.
Reactive parameters can also be modified:
--r_eq_12 R_EQ_12 r_eq for reactive atoms (1-2 interaction)
--eps_12 EPS_12 epsilon for reactive atoms (1-2 interaction)
--r_eq_13_scaling R_EQ_13_SCALING
r_eq scaling for 1-3 interaction across reactive atoms
--r_eq_14_scaling R_EQ_14_SCALING
r_eq scaling for 1-4 interaction across reactive atomsReceptor preparation can't handle heteroatoms for the moment. Also nucleic acids, ions, and post-translational modifications (e.g. phosphorilation) are not supported. Only the 20 standard amino acids can be parsed, and it is required to have Amber atom names and hydrogens. No atoms can be missing.
Make affinity maps for the _rigid.pdbqt part of the receptor.
Make affinity maps for the _rigid.pdbqt part of the receptor. mk_prepare_receptor.py will prepare the GPF for you.
mk_prepare_ligand.py -i sufex1.sdf --reactive_smarts "S(=O)(=O)F" --reactive_smarts_idx 1 -o sufex1.pdbqt\
For reactive docking there are two options that need to be passed to AutoDock-GPU:
For reactive docking there are an additional option that needs to be passed to AutoDock-GPU:
console --import_dpf
The --derivtype option, if needed, was written by mk_prepare_receptor.py to a file suffixed with .derivtype.
The filename to be passed to --import_dpf was written by mk_prepare_receptor.py
and it is suffixed with reactive_config.
ADGPU -I *.reactive_config -L sufex1.pdbqt -N sufex1_docked_ -F *_flex.pdbqt -C 1