Use residue template matching to prepare ligands in PDB format

[diogomart]

Kudos to @rwxayheee for this idea. Would enable:

mk_prepare_ligand.py -i ligand.pdb

as long as the residue name can be looked up in the CCD

[rwxayheee]

Hi @diogomart

I can work on this some time later. Just a few thoughts:

It's easy to do with Python scripting (see example at the end), but adding it to mk_prepare_ligand.py requires a PDB parser or Polymer constructor, and ResidueChemTemplates. Currently, the CLI script is designed more for batch processing of ligands from SDF, not so optimal for covalent ligands or ligands from PDB. Part of it can be re-used, while there are other things that are different from the standard ligand preparation.

It needs a little more work and the addition could make the CLI script even longer. I would also like to add options to modify/override a template from CCD.

Here's a Python example for people interested in trying. No input is required. Prody fetches the PDB file with the given PDB ID, and the ligand is specified by residue name.

#!/usr/bin/env python

""""
micromamba create -c conda-forge -n meeko_tutorial_py39 python=3.9 -y
micromamba activate meeko_tutorial_py39

micromamba install -c conda-forge numpy scipy rdkit gemmi meeko -y

pip install prody
"""

from rdkit import Chem

import prody

from meeko import MoleculePreparation
mk_prep = MoleculePreparation()
from meeko import PDBQTWriterLegacy
from meeko import Polymer
from meeko import ResidueChemTemplates

# PDB ID or PDB filename of the starting ligand structure
pdb_token = "4kqy"

# ProDy selection language to define ligand atoms
lig_atoms_selection = "resname SAM"

pdb_prody_mol = prody.parsePDB(pdb_token)
lig_prody_mol = pdb_prody_mol.select(lig_atoms_selection)

templates = ResidueChemTemplates.create_from_defaults()
polymer = Polymer.from_prody(
    lig_prody_mol,
    templates,
    mk_prep,
    default_altloc="A",
)

output_bname = "lig_"

for lig_resid in polymer.monomers.keys(): 
    resid = lig_resid
    mol = polymer.monomers[lig_resid].rdkit_mol # should have all Hs
    atom_names = polymer.monomers[lig_resid].atom_names
    for i,atom in enumerate(mol.GetAtoms()): 
        pdb_residue_info = Chem.AtomPDBResidueInfo(atomName=atom_names[i])
        atom.SetMonomerInfo(pdb_residue_info)

    molsetup_list = mk_prep(mol)
    molsetup = molsetup_list[0]

    pdbqt_string = PDBQTWriterLegacy.write_string(molsetup)[0]

    output_pdbqt = output_bname + f"{pdb_token}_{resid}.pdbqt".replace(":","-")
    with open(output_pdbqt, 'w') as f:
        f.write(pdbqt_string)
    print(f"File written: {output_pdbqt}")

[diogomart]

Either a new script or adding to mk_prepare_ligand.py look good to me 👍