PDB reading problem

[avnikonenko]

Dear authors, I tried to read with RDKit pdbqt file generated by the latest dev version of Meeko, but I get an error about 27 atom of H which as I understand is because it is shifted 1 position to the left in the columns. Also Meeko for some reason wrote SER residue name instead of UNK for some atoms in this generated pdbqt file. Files: error_H.mol.txt error_H.pdbqt.txt Code:

>> import rdkit
>> rdkit.__version__
'2021.03.5'
>> from rdkit import Chem
>> mol = Chem.MolFromMolFile('error_H.mol.txt', removeHs=False)
>> mol
<rdkit.Chem.rdchem.Mol object at 0x7fa1206d74e0>
>> from openbabel import openbabel as ob
>> from meeko import MoleculePreparation
>> from meeko import obutils
>> mol_ob = obutils.load_molecule_from_string(Chem.MolToMolBlock(mol), molecule_format='MOL')
>> preparator = MoleculePreparation(merge_hydrogens=True, macrocycle=False, hydrate=False, amide_rigid=True)
>> preparator.prepare(mol_ob)
>> pdbqt_string = preparator.write_pdbqt_string()
>> mol_pdbqt = Chem.MolFromPDBBlock('\n'.join([i[:66] for i in pdbqt_string.split('\n')]))
[15:53:42] Cannot determine element for PDB atom #27

Could you help with it please?

[diogomart]

The SER residues are assigned by OpenBabel. It seems that hydroxyls sitting two carbons away from an amide are at risk of being classified as SER. OpenBabel also adds some extra spaces to some of the hydrogen atom names:

>>> import openbabel
>>> print(openbabel.__version__)
3.1.0
>>> from openbabel import openbabel as ob
>>> obconv = ob.OBConversion()
>>> obconv.SetInFormat('mol')
True
>>> obmol = ob.OBMol()
>>> obconv.ReadFile(obmol, 'error_H.mol.txt')
True
>>> for atom in ob.OBMolAtomIter(obmol):
...     res = atom.GetResidue()
...     print('atom_name: "%s", res_name: "%s", res_num: %d' % (res.GetAtomID(atom), res.GetName(), res.GetNum()))
... 
atom_name: "N ", res_name: "UNK", res_num: 0
atom_name: "C ", res_name: "UNK", res_num: 0
[...]
atom_name: "H ", res_name: "UNK", res_num: 0
atom_name: "H  ", res_name: "SER", res_num: 1
atom_name: "HA ", res_name: "SER", res_num: 1
atom_name: "HB1", res_name: "SER", res_num: 1
atom_name: "HB2", res_name: "SER", res_num: 1
atom_name: "H ", res_name: "SER", res_num: 1
atom_name: "H1", res_name: "SER", res_num: 1
atom_name: "H2", res_name: "SER", res_num: 1
atom_name: "H  ", res_name: "SER", res_num: 1
atom_name: "H ", res_name: "UNK", res_num: 0
[...]

In any case, a better alternative to building RDKit Mol objects from PDB files (which is dangerous because all bonds will be single bonds and then benzenes become cyclohexanes) is to do the following:

import rdkit
from rdkit import Chem
print(rdkit.__version__)

mol = Chem.MolFromMolFile('error_H.mol.txt', removeHs=False)

from openbabel import openbabel as ob
from meeko import MoleculePreparation
from meeko import obutils
mol_ob = obutils.load_molecule_from_string(Chem.MolToMolBlock(mol), molecule_format='MOL')
preparator = MoleculePreparation(merge_hydrogens=True, macrocycle=False, hydrate=False, amide_rigid=True)
preparator.prepare(mol_ob)
pdbqt_string = preparator.write_pdbqt_string()

obconv = ob.OBConversion()
obconv.SetOutFormat('mol')

tmpfn = '.temporary-file-to-init-PDBQTMolecule.pdbqt'
with open(tmpfn, 'w') as f:
    f.write(pdbqt_string)

from meeko import PDBQTMolecule
pdbqt_mol = PDBQTMolecule(tmpfn, is_dlg=False) # forgot to implement init from string, sorry!
mols_pdbqt = []
for pose in pdbqt_mol: # just 1 pose here, but this will work for docking outputs
    copy_obmol = ob.OBMol(mol_ob) # need a fresh copy for each new pose
    pose.copy_coordinates_to_obmol(copy_obmol)
    mol_pdbqt = Chem.MolFromMolBlock(obconv.WriteString(copy_obmol))
    mols_pdbqt.append(mol_pdbqt)

[avnikonenko]

Thank you so much!