Adding Hydrogen without changing the ligand pose

[Furaxixi]

Hi Molsrub Team,

Thank you for your amazing work.

I have a ligand in an SDF file, and I’d like to use scrub.py to add hydrogens at a specific pH without altering the existing conformation/pose.

I’ve tried different commands, including --skip_gen3d, but it either generates a 2D representation of the ligand or alters the conformation.

Could you advise on the correct approach to achieve this?

Looking forward to your response.

Kind regards,
Francis

[rwxayheee]

Hi @Furaxixi

RDKit's AddHs could be the most basic and convenient option if you only want to add hydrogens based on the protonation states computed from molscrub. But there are some limitations of AddHs to consider, including not having the hydrogens added at the optimal places and not accommodating the necessary conformation changes that are required by the transitions of protonation or tautomeric states, etc. If you're interested, see the example at the end for a single amino acid Glutamine. In the outputs, the carboxylic acid group is protonated in the anti-configuration.

Therefore, this simple approach does the minimal task but is not 100% correct. Depending on what kind of analysis you're proceeding with the output structures, you may find the limitations acceptable, or might need to do some further structural optimization. These output structures should be be ok as starting points for geometry optimization or dynamics simulation, but may not be suitable for energy calculations that doesn't include optimization of the hydrogens.

Get the initial conformer from:
https://pubchem.ncbi.nlm.nih.gov/compound/5961#section=3D-Conformer

from rdkit import Chem
from rdkit.Chem import AllChem
from scrubber import Scrub
from pathlib import Path

# Molecule from SDF
sdf_fn = "Conformer3D_COMPOUND_CID_5961.sdf"
sdf_mol = Chem.MolFromMolFile(sdf_fn, removeHs=True)
smi = Chem.MolToSmiles(sdf_mol)

# Molecule input for Scrubber
mol = Chem.MolFromSmiles(smi)
# Scrubber configuration
scrub = Scrub(
    ph_low=4.5,
    ph_high=10.5,
    skip_gen3d=True,
)

# Enumerate output protonation states
for i, mol_state in enumerate(scrub(mol)):
    
    # Generate the template molecule with scrubbed protonation state
    srubbed_smi = Chem.MolToSmiles(mol_state, isomericSmiles=True, canonical=True)
    template_mol = Chem.MolFromSmiles(srubbed_smi)
    Chem.SanitizeMol(template_mol)
    
    # Use the template molecule to assign bond orders in molecule copy from SDF
    new_mol = AllChem.AssignBondOrdersFromTemplate(template_mol, sdf_mol)
    new_mol = Chem.AddHs(new_mol, addCoords=True)
    
    # Write the new molecule to SDF
    new_fn = f"{Path(sdf_fn).stem}_{i}.sdf"
    writer = Chem.SDWriter(new_fn)
    writer.write(new_mol)
    writer.close()
    print(srubbed_smi, "->", new_fn)