mda equivalent of Rdkit's addHs

[corey-taylor2]

G'day all,

Looking to extract a ligand and surrounding residues from a Universe but then to cap valences in a manner similar to Rdkit's addHs module to feed into some qchem workflows. Directly saving the selection as a pdb shows there are indeed unsatisfied valences.

I've Googled and searched and can't find a way to do this, so either it's so blindingly obvious that this topic will prove a useful reminder to others.... Or it's not as straightforward as I'd like it to be. Anyone have some code laying around to achieve my dreams? Cheers all.

[hmacdope]

You can use the RDKit converter to go directly to RDKit and then add the Hs there!

from rdkit import Chem
import MDAnalysis as mda
u = mda.Universe("protein.pdb")
atoms = u.select_atoms("name CA")
rdkit_mol = atoms.convert_to("RDKIT")
# do whatever you like in RDKit.

[hmacdope]

Yep should work!

[orbeckst]

@corey-taylor2 if this works for you, can you mark @hmacdope 's answer as this will help other readers, too (and I am sure it will make @hmacdope happy)? If it doesn't work, let us know why, that's also valuable feedback!

[corey-taylor2]

Yup, done!

[corey-taylor2]

Update: there might be a missing step or two with the conversion?

import MDAnalysis as mda
import psi4
import rdkit.Chem

example_complex = 'complex.pdb'
complex = mda.Universe(example_complex)
protein_selection = complex.select_atoms("resname LIG or ((protein) and byres around 6.0 group ligand)", ligand=ligand)
rdkit_mol = protein_selection.convert_to("RDKIT")

Doing this, I get a NameError:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[27], line 1
----> 1 rdkit_mol = protein_selection.convert_to("RDKIT")

File ~/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/core/accessors.py:204, in ConverterWrapper.__call__(self, package, *args, **kwargs)
    [201](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/core/accessors.py?line=200) except AttributeError:
    [202](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/core/accessors.py?line=201)     raise ValueError(f"No {package!r} converter found. Available: "
    [203](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/core/accessors.py?line=202)                      f"{' '.join(self._CONVERTERS.keys())}") from None
--> [204](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/core/accessors.py?line=203) return convert(*args, **kwargs)

File ~/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py:359, in RDKitConverter.convert(self, obj, cache, NoImplicit, max_iter, force)
    [357](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=356) kwargs = dict(NoImplicit=NoImplicit, max_iter=max_iter, force=force)
    [358](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=357) if cache:
--> [359](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=358)     mol = atomgroup_to_mol(ag, **kwargs)
    [360](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=359)     mol = copy.deepcopy(mol)
    [361](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=360) else:

File ~/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py:428, in atomgroup_to_mol(ag, NoImplicit, max_iter, force)
    [426](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=425) # attributes accepted in PDBResidueInfo object
    [427](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=426) pdb_attrs = {}
--> [428](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=427) for attr in RDATTRIBUTES.keys():
    [429](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=428)     if hasattr(ag, attr):
    [430](file:///home/ubuntu/miniconda3/envs/qchem/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py?line=429)         pdb_attrs[attr] = getattr(ag, attr)

NameError: name 'RDATTRIBUTES' is not defined

What am I missing in my inputs? My s/w versions, if that's relevant.

rdkit             2024.9.4
MDAnalysis        2.8.0

[corey-taylor2]

If it helps, I get the same error if I run the following toy example:

mda.Universe.from_smiles("CCO").atoms.convert_to("RDKIT")